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SUMMARY 


This study concerns the rotordynamic characteristics of fluid-encompassed rotors, 
with special emphasis on shrouded pump impellers. The core of the study is a versa- 
tile and categorically new finite-element-based perturbation model, which is based on a 
rigorous flow analysis and what we have generically termed the “virtually” deformable 
finite-element approach. The model is first applied to the case of a smooth annular seal 
for verification purposes. The rotor excitation components, in this sample problem, give 
rise a purely cylindrical, purely conical and a simultaneous cylindrical/conical rotor whirl 
around the housing centerline. In all cases, the computed results are compared to ex- 
isting experimental and analytical data involving the same seal geometry and operating 
conditions. Next, two labyrinth-seal configurations, which share the same tooth-to-tooth 
chamber geometry but differ in the total number of chambers, were investigated. The 
results, in this case, are compared to experimental measurements, for both seal configu- 
rations. The focus is finally shifted to the shrouded-impeller problem, where the stability 
effects of the leakage flow in the shroud-to-housing secondary passage are investigated. To 
this end, the computational model is applied to a typical shrouded-impeller pump stage, 
fabricated and rotordynamically tested by Sulzer Bros., and the results compared to to 
those of a simplified “bulk-flow” analysis and Sulzer Bros.’ test data. In addition to as- 
sessing the computed rotordynamic coefficients, the shrouded-impeller study also covers a 
controversial topic; namely that of the leakage-passage inlet swirl, which w’as previously 
cited as the origin of highly unconventional (resonance-like) trends of the fluid-exerted 
forces. In order to validate this claim, a “microscopic” study of the fluid/shroud interac- 
tion mechanism is conducted, with the focus being on the structure of the perturbed flow 
field associated with the impeller whirl. The conclusions of this study were solidified by 
the outcome of a numerical-certainty exercise, where the grid dependency of the numeri- 
cal results is objectively examined. The final phase of the shrouded-impeller investigation 
involves the validation of a built-in assumption, in all other perturbation models, whereby 
single-harmonic tangential distributions of all the flow thermophysical properties are im- 
posed. Grid dependency of the fluid-induced The last phase of the investigation course is 
aimed at verifying the fine details of the perturbed flow field in light of recent set of de- 
tailed LDA measurements in a smooth annular seal. Grid dependency of the fluid-induced 
forces is also investigated, and specific recommendations made. 



OVERVIEW AND REPORT STRUCTURE 


This report presents the development, verification, and applications of a categori- 
cally unprecedented perturbation approach to the rotordynamic problem arising from the 
fluid/rotor interaction in turbomachinery applications. This is a finite-element-based ap- 
proach, which is based on the perception that the fluid reaction forces developed in the 
rotor-to-housing gap as a result of the rotor whirl, arise from infinitesimally small defor- 
mations of the finite elements in this gap. Although the research tasks were all focused on 
the rotordynamic characteristics of shrouded impellers (Fig. la), the problem formulation 
is sufficiently general to embrace other conceptually similar problems, all belonging to the 
isolated-seal category. 

Since many aspects of the current study have already been published, or submitted 
for publication, and for the sake of compactness, the reader is referred to the manuscript 
of a specific paper, in one of the two appendices, each time reference is made to the subject 
of this paper. These two appendices contain already published journal articles, and those 
which are currently under review, respectively. Once the reader is referred to one of these 
papers, only the relevant contents of the paper are recited in the core of the report, should 
that serve the purpose of clarity. The appendices, in this sense, are hardly peripheral, but 
constitute an equally essential segments of the text. 

Among the segments which were overridden in the report core is that concerning 
the mathematical development of the “virtually” deformable finite-element-based pertur- 
bation model. This tedious and predominantly mathematical subject is presented in the 
subsection “Paper# 1” of Appendix 1 for the reader who might be interested in the mathe- 
matical foundation of the model, and the finite-element formulation of the general problem 
of fluid/rotor interaction. Other papers of the appendices involve exhaustive discussions 
of certain applications of the model, and as such, may distract the user from the main 
application focus, being that of shrouded impellers, had they been included in the report 
core. Of these, Paper# 2 concerns the model application to the case of a smooth annu- 
lar seal with a cylindrically whirling rotor, and confirms the traditional assumption of a 
single-harmonic variation of the shroud pressure perturbation around the circumference. 
This assumption is indicative of a set of restrictions that has been (almost blindly) common 
among all of the existing models in this engineering discipline, regardless of the physics 
of the problem (the reader is alerted to the implied fact here, that such an assumption is 
not made in the current computational model). Next, Paper# 3, in the same appendix, 
expands the model to the unstable seal operation mode involving a conical rotor whirl, 
again in a smooth annular seal. The two degrees-of-freedom of the rotor axis, meaning 
the linear and angular displacements, are then lumped together and applied to the same 
seal configuration. The resulting cylindrical and conical whirl components of the rotor 
axis were analyzed as occurring simultaneously, and the results appraised in Paper# 4 of 
Appendix 1. The last paper (# 5) in this appendix marks the beginning of the process 
where the unstable operation of a shrouded pump impeller is investigated. This paper 
presents the highly irregular features of the unperturbed (zeroth-order) flow fields in two 
typical leakage-passage configurations, with two, conventional wear-ring and face, seals 
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comprising the passage tight-clearance regions. A comprehensive comparison between the 
leakage (or secondary) flow fields, associated with these seal categories, is made in this 
paper, with a special emphasis on rotordynamics-related aspects such as the leakage- 
passage inlet swirl angle. 

Appendix 2 in this report contains three recent papers which are currently at different 
phases of the review process for publication. The first of these papers is a natural con- 
tinuation of Paper# 4 in Appendix 1. This paper contains the rotordynamic coefficients 
associated with one of the two seal configurations mentioned above, namely the face-seal 
configuration. The paper also contains a comparison between the computed coefficients 
versus a bulk-flow solution (Childs, 1989) and the experimental data of Sulzer Bros. (Bol- 
leter et al., 1989). Next, the current perturbation approach is applied to the category of 
labyrinth seals, where the basic issue of the tooth-to-tooth chamber count impact on the 
seal stability is investigated. The last paper (# 3) in Appendix 2 makes use of recent LDA 
measurements, by Morrison et al. (1992), in assessing the computed perturbations of the 
flow-velocity components throughout an annular seal with a cylindrically whirling rotor. 
The paper shows that there indeed is a notable agreement between the two sets of results 
in the downstream seal segment, and that the most observable agreement is progressively 
attained as the seal-discharge station is approached. This, as stated in the paper, is a com- 
prehendable finding, in the sense that this downstream segment of the seal is where the real 
admission effects, in the actual rig, have diminished, allowing a one-to-one comparison of 
the computed and measured sets of data. 

The core segment of this report addresses a certain critical issue, all related to the 
major application of the current computational model, namely that of the shrouded pump 
impeller. The issue here was raised by Childs (1989), who that an inlet swirl ratio (leakage- 
passage inlet swirl component/impeller tip speed) that is above 0.5 is sure to produce 
abrupt fluctuations in the fluid reaction forces in a narrow range of positive whirl ratios. 
Childs’ conclusions were based on a simplified bulk-flow analysis, and perhaps are limited 
to the Sulzer Bros.’ pump impeller. Since such fluctuations would greatly influence the 
impeller rotordynamic coefficients, as well as prohibit the use of a simple quadratic inter- 
polation of the fluid-exerted force components, at least in the disputed whirl frequency 
range, a thorough investigation of this phenomenon was warranted. Lack of conclusiveness, 
meanwhile, was caused by Guinzburg’s experimental results (1992) which dismissed Childs’ 
observation. Although Guinzburg’s rig design did not allow, or even simulate, any primary 
flow effects, the categorical denial of this “resonance-like” phenomenon was, at least, sur- 
prising. Seriousness of this issue was further established, when a careful review of Sulzer 
Bros.’ experimental measurements revealed that the fluid-forces versus whirl-frequency 
indeed exhibit some fluctuations (Fig. lb), although the latter were much milder by com- 
parison. Vagueness of the entire issue was later concluded, as the preliminary results of the 
current perturbation model came to support the existence of such force fluctuations, with 
the exception that the severity of these fluctuations was closer to those of Sulzer Bros, 
than they were to Childs’ force trends. A full-scale investigation of the issue was then 
conducted. It is exclusively this investigation which is emphasized in the core segment of 
this report. 


Investigation of the fluid-exerted forces in the Sulzer Bros.’ pump case involved 
several potentially-contributing factors concerning the current model as well as the bulk- 
flow model. These factors were as follows: 

1. Grid dependency of the fluid-exerted forces in our analysis (Fig. 2). 

2. The manner in which the shroud pressure perturbations are altered during the transition 
from negative (backward) to positive (forward) whirl. 

3. Effect of the leakage passage inlet swirl velocity on our computed forces. 

4. Validity of the bulk-flow-model assumption that the perturbations of the flow properties 
(particularly the shroud pressure) are expressible via a single harmonic in the circumfer- 
ential direction. 

The last section of the report-core segment concerns proposed research tasks which 
are aimed at continuing the code validation process, on one hand, and expanding the 
current perturbation model towards either more prediction accuracy or new applications, 
on the other. One of the proposed topics is a bench-mark test case, and involves numerical 
simulation of the Cal. Tech, test rig (Guinzburg, 1992). This rig was designed with the 
intention of testing the rotordynamic consequences of the isolated leakage passage. The 
passage inlet swirl in this case was varied from -200% to +200% of the rotor tip speed, and 
no indication of force fluctuation was observed throughout the entire inlet swirl range. 

Also proposed is another related study, which has to do with the flow domain with 
which the shroud interacts. As seen in Fig. la, we have so far simulated the pri- 
mary/leakage flow interaction at the impeller inlet and discharge stations. This is ob- 
viously an upgrade to the common approach where the leakage passage is totally isolated, 
and fictitiously uniform boundary conditions are specified at both ends. Nevertheless, we 
believe that a much more realistic model of the fluid/shroud interaction can be attained 
by physically including the impeller primary passage as part of the computational domain. 
In this case, contributions to the shroud forces will (as they should) be produced in both 
flow passages. Implementation of this upgrade will undoubtedly improve the accuracy of 
the computed rotordynamic coefficients. 

A major expansion of our perturbation model is that concerning the hydrostatic bear- 
ing rotordynamics problem. Numerical models which were developed over the past three 
years were generally based on the bulk-flow theory in which inertia domination (at the 
recess edges and in the film land) were simulated through simple one— dimensional means, 
and the flow turbulence was recognized only at the walls. These models, more or less, 
uncoupled the recess and film-land contributions to the journal forces and, consequently, 
treated the fluid/journal interaction problem in a steady-state sense. Our approach, on 
the other hand, is that of repeatedly marching in time through a circumferential “pitch” 
between two successive recesses, and computing the instantaneous fluid-exerted forces dur- 
ing the process until until convergence (in terms of a fixed cyclic force pattern) is achieved. 
Moreover, contrary to the existing models, where no journal-to-housing flow gradients 
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Fig. 2 Dependency of the fluid-exerted forces on the grid resolution in the circumferential direction 
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are taken into account, the proposed model is based on the fully three-dimensional flow- 
governing equations, coupled with a provenly accurate combination of turbulence closure 
and near-wall flow analysis. 
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COMPUTATIONAL PROCEDURE 


Analysis of the fluid/shroud interaction within a pump stage, is based on what we in- 
terimly termed the “virtually” deformable finite— element concept (Baskharone and Hensel, 
1991). Under this approach, perturbations in the fluid thermophysical properties including, 
in particular, the shroud pressure are perceived as a direct result of varied distortions in the 
finite element assembly occupying the shroud-to-housing gap as the impeller undergoes a 
whirling motion. 

The computational process is initialized by solving the centered-impeller flow field 
throughout the computational domain shown in Fig. 3. This “zeroth-order” flow solution 
is then used as input to the perturbation model. This step physically corresponds to the 
point when the impeller whirl (defined by a virtual impeller eccentricity e and a finite whirl 
frequency $7) is recognized (Fig. 4). This is also when the computational domain becomes 
totally three-dimensional due to the impeller eccentricity. The finite— element discretization 
unit, in this case, is the twenty-noded curve— sided isoparametric element shown in Fig. 5, 
with the frame of reference being attached to the impeller and whirling with it (Fig. 6) 
in order to eliminate the time-dependency of the shroud-to— housing flow field. Expansion 
of the finite-element equations (on an element-by-element basis) in terms of the nodal 
values of displacement, eventually gives rise (once assembled into the global system) to an 
equivalent set of equations the unknowns in which are [dPi/de, i = 1,2,..., A], where pi 
is a general flow property (e.g velocity component or pressure) at the computational node 
“i” . Of these, the shroud pressure perturbations are extracted and integrated over the 
entire shroud surface, giving rise to the rate (with respect to e) at which the fluid forces 
on the shroud are developed. These forces are, in the end, used to compute the direct and 
cross-coupled stiffness, damping and inertia (added mass) coefficients of the fluid/shroud 
interaction. 
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Quadrilateral curve-sided isoparametric finite element for analyzing the perturbed flow field 







UNIQUENESS AND IMPACT OF THE NEW APPROACH 

Judging by the initial impression of a typical discusser (e.g. a paper reviewer, a 
conference attendee,. ..etc), the very foundation of the current perturbation model can 
easily be misconceived. For instance, referring to the use of the finite-element method (in 
conjunction with the model) as a means of solving the “perturbed flow field” or solving 
“the perturbation equations” is misleading and totally wrong, respectively. The fact (by 
reference to Paper# 1 in Appendix 1) is that the perturbed-field equations themselves 
arise from physical displacements, velocities and acceleration components which are viewd 
from a moving (rotating and translating) frame of reference, as a result of the rotor whirl 
(Fig. 6). These are carefully superimposed on the unperturbed flow field, giving rise 
to a new system (governing equations and boundary conditions) under which the fluid 
response is totally unrestricted . When this approach is compared to existing perturbation 
models (where the circumferential variation of the fluid response is externally imposed as 
a single harmonic), the question becomes: “what if the flow physics do not naturally lend 
themselves to such pre-specifed response”, or, “what if the circumferential variation of the 
clearance width (itself) is not a single harmonic^ such as the case of a tagentially-deformed 
shroud or a hydrostatic bearing)”. 

The process of “adding” whirl-related terms to the unperturbed flow field (above) 
involves all of the flow kinematical properties, with the common factor being the rotor 
eccentricity “e”. For instance, the added displacements are those of the nodal points 
(relative to an observer at the origin of the whirling frame of reference, Fig. 4), the added 
velocity components are, for instance, those of a node on the housing surface (since the 
housing is non-stationary in this frame of reference, and the added acceleration components 
are those of the Coriolis and Centriputal types, and appear as a result of the origin motion, 
since the latter is attached to the rotor (Fig. 4). 

The problem formulation (above) is hardly a mere mathematical alteration of existing 
perturbation models. The critical issue here is “whether one should define the rotor eccen- 
tric motion and the fluid response to it.” To the author, the answer is “not always, unless 
the problem is the simplest possible (e.g. the case of a straight annular seal).” In fact, 
such a pre-imposed response (namely that of the single-harmonic type) is proven (later in 
this report) to produce an error that is as high as 21% within regions in the leakage path 
of a pump, where such an assumption is inconsistent with the local flow physics. 

In view of the above, a traditional question calling for a comparison between the 
finite-element technique and other methods, is perhaps less relevant to this study than the 
asker might think. 
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RESULTS AND DISCUSSION 


Origin of the Rotordynamic Force Fluctuations: 

The study outlined in Appendix 1, Paper# 4 concerns the centered-impeller flow field 
in the Sulzer Bros, face— seal pump configuration, and compares it to a twin pump geometry 
where the leakage-control device is of the conventional wear-ring type. Of particular 
interest in this paper is the swirl velocity distribution in the secondary— flow passage leading 
to the seal (Figs. 6 & 9 of the Appendix). Comparison of the seal-inlet swirl velocity ratios, 
in both pump configurations, with those of Childs (1989), revealed one possible reason why 
the shroud forces in Childs’ study had the nearly resonance-type trends. The reason, we 
believe, is that the seal— inlet swirl velocity (obtained through simplified matching means 
in Childs’ analysis) was approximately double our computed values. Note that this swirl 
ratio in our analysis is an internal variable that is totally unrestricted. 

Fig. 2 contains our rotordynamic results for the face-seal pump configuration using 
the centered-impeller flow solution, mentioned above, as input. The results in this Fig- 
ure was particularly revealing to us. The reason is that initially— obtained shroud forces 
were the outcome of a finite— element model with only seven computational planes in the 
circumferential direction. This, by reference to the Figure, is much less than the number 
required for grid independency. 


Effect of the Leakage-passage Inlet Swirl 

A test case, involving a modified version of Sulzer Bros.’ face-seal pump configuration 
(Fig. 3), was investigated. The intention here was to investigate the effect of the leakage- 
passage inlet swirl on the impeller rotordynamic characteristics. In this case, the inlet swirl 
was totally destryed, a situation that is physically achievable by inserting a swirl brake at 
the leakage-passage inlet station, or through the use of a cascade of deswirl vanes in the 
primary passage downstream from the impeller (Fig. 7). 

Results of the centered-impeller flow analysis, in the absence of leakage-passage inlet 
swirl, are shown in Figs. 8,9 &10. These are vector and contour plots of the velocity 
and pressure fields. The plots in this Figure were then compared to their counterparts 
in Appendix 1, Paper# 4, with the latter corresponding to the design-point value of the 
leakage-passage inlet swirl (Figs. 8,9 &10 in Appendix 1, Paper# 4). The conclusion 
was that “killing” the inlet swirl had no appreciable effect on the meridional velocity 
vectors throughout the leakage passage. This was, to some extent, anticipated since the 
recirculatory motion in this passage is primarily controlled by two other factors; namely the 
fluid migration radially outwards (near the shroud) as a result of the centrifugal force, and 
inwards (near the housing) due to the streamwise static pressure differential. However, the 
swirl-velocity/tip-speed ratio shown in Fig. 9 seems to remain much smaller throughout 
the leakage passage by comparison. This, to us, was the most significant difference between 
the two sets of results, and was regarded as a sufficient reason to expect significant 
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Fig. 7 Deswirl vanes as a means of eliminating the leakage-passage inlet swirl 
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Fig. 9 Contours of the nondimensional swirl velocity for the zero leakage-passage preswirl case 
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improvements in the stability-related rotordynamic coefficients. As for the static pressure 
distributions (Fig. 10 and Fig. 10 in Appendix 1, Paper# 4), no major gain seems to have 
been obtained by destroying the leakage-passage inlet swirl velocity, with the exception 
that the secondary passage segment leading to the seal now shares the streamwise pressure 
differential with the seal, more so than the case was under the design-point preswirl (Fig. 

10 in Appendix 1, Paper# 4). This feature is desirable from a design standpoint in the 
sense that utilization of a lower-resistance seal would, in this case, be allowed. Finally, 
the swirl velocity profiles slightly downstream from the leakage-passage inlet station (Fig. 

11 and Fig. 5 in Appendix 1, Paper# 4) were compared. As expected, a rather thick 
boundary-layer type region exists near the shroud, under the no-preswirl pump operation, 
where the swirl velocity is notably high as a result of excessive flow recirculation in this 
region. Note that the swirl velocity component in Figs. 9 & 11 is non-dimensionalized 
using the impeller tip speed. The non-dimensional pressure in Fig. 10, on the other hand, 
is defined as follows: 

_ ( p - Pi 

P= r 2 

Pit 

where p and p t are the local and stage-inlet static pressures, respectively, p is the fluid 
density and Ut is the impeller tip speed. 

The “no-preswirl” zeroth-order flow solution was then used to compute the fluid- 
exerted forces over a range of whirl frequencies between -125% and +125% of the impeller 
speed. The results are shown in Fig. 12, together with the computed forces and experi- 
mental data, both of which corresponding to the design-point value of the leakage-passage 
inlet swirl. The first impression one would have by examining the tangential-force curve 
in this Figure is that the no-preswirl pump operation is comparatively more stable (as 
anticipated) judging by the smaller range of positive whirl frequency within which the 
tangential force is positive. 

The impeller rotordynamic coefficients, corresponding to zero leakage-passage preswirl, 
are compared to the numerical and experimental data, associated with the design-point 
preswirl, in Table 1. All three sets of data in this table were obtained through a quadratic 
least-square fit of the fluid-exerted forces, as outlined by Baskharone and Hensel, 1991. 
Referring to Table 1, and limiting the comparison to the computed coefficients (last two 
columns of the Table), it is clear that the zero leakage-passage preswirl gives rise to a 
much stable impeller operation. While the direct damping coefficient “C” remains almost 
unaffected, the cross-coupled stiffness coefficient “k” becomes negative (thus stabilizing) 
in the zero-preswirl operation mode. 


Examination of the Perturbed Flow Field: 

In an attempt to understand the origin of rotordynamic forces, we developed a post- 
processor which, in effect, “unwraps” the shroud surface, in the manner illustrated in 
Fig. 13, and displays the distribution of pressure perturbations (i.e. dp/de over the 
“unwrapped” surface, where “e” is the impeller eccentricity. Note that integration of the 
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ig. 11 Profile of the nondimensional swirl velocity near the leakage— passage inlet station 
for the zero-preswirl case 


EXPERIMENTAL DATA (WITHOUT DESWIRL VANES) 
CURRENT ANALYSIS (WITHOUT DESWIRL VANES) 
CURRENT ANALYSIS (WITH DESWIRL VANES) 



Fig. 12 Comparison of the fluid exerted forces with and without deswirl vanes 
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pressure perturbations over the entire shroud surface yields the net fluid-exerted force 
which, in turn, is resolved in the radial and tangential directions. The pump operation 
with and without leakage-passage preswirl were both under consideration in a range of 
whirl frequency ratios between -0.2 and +0.3. The main objective here was to see whether 
there are any redistribution of the pressure perturbation contours as the impeller shifts 
from slightly negative to slightly positive whirl frequencies. 

Fig. 14 (on three successive pages) shows the pressure perturbation contours in the 
selected whirl-frequency range. For backward whirl (first two whirl frequency ratios), 
the Figure shows that the zero-presw T irl pressure perturbations peak, more or less, at the 
minimum-clearance position (circumferential position c-c in Fig. 13). The Figure also 
shows that the design-point preswirl pump operation gives rise to pressure perturbation 
peaks which are clearly ahead of the minimum-clearance position in the direction of rotation 
The relatively symmetric distribution of the pressure perturbation in the zero-preswirl case 
is perhaps the reason why the tangential force in Fig. 12 is comparatively small in this 
negative whirl-frequency range. 

Proceeding to the positive whirl-frequency ratios in Fig. 14, it is seen that both pump 
operation modes give rise to peak pressure perturbations ahead of the minimum-clearance 
position in the direction of rotation. This occurs in the whirl frequency ratio range between 
0.0 and +0.2, and indeed represents a shift in the case of the zero-preswirl case. This shift 
causes the notably rise in the tangential force curve (Fig. 12) in this frequency range. 
While the pressure perturbation peak in Fig. 14 continues to exist ahead of the minimum- 
clearance position (c-c) for the pump operation with inlet preswirl, a rather surprising 
contour redistribution occurs for the zero-preswirl operation, whereby the pressure per- 
turbation peak is now at a circumferential location that is lagging the minimum-clearance 
position (relative to the direction of rotation) at fl/u> = +0.3. Referring to Fig. 12, note 
that this shift causes the tangential force (for zero-preswirl) to be negative and, therefore, 
stabilizing. 

The foregoing discussion illustrates the significance of the peak pressure-perturbation 
position (on the shroud surface), but lacks experimental validation. Fortunately, a recent 
set of LDA velocity measurements by Morrison et al. (1992) made a qualitative verification 
of our results possible. Fig. 15 shows a typical contour plot of the non-dimensional 
through-flow velocity which was obtained in this experimental study for an annular seal 
with a synchronously whirling rotor (fl/u; = 1.0. Shown in the same Figure is our computed 
non-dimensional axial-velocity perturbation contours in the leakage passage at the location 
marked x-x for a synchronous whirl-frequency ratio. Note that this Figure does not allow 
a one-to-one comparison, for the simple reason that Morrison’s contours are those of the 
non-dimensional through-flow velocity: 

V z 

v Ziin 

while our contours are those of the non-dimensional perturbation: 

dV z 
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where V z> i n is the seal-inlet through-flow velocity (in Morrison’s rig), and c is the local 
clearance at the location x-x (in our case). Nevertheless, The Figure provides a means for 
qualitative comparison, particularly in regards to the tangential position of the peak value 
of through-flow velocity. Both the experimental and computed contours in Fig. 15 clearly 
indicate that the maximum disturbance position (marked “H” in Fig. 15) is at a tangential 
location in the opposite direction (relative to the rotor speed) from the minimum-clearance 
position. Pursuing this qualitative comparison, the shroud-surface pressure perturbation 
contours were also generated, and are shown in Fig. 16 (also for a synchronously-whirling 
impeller). Examination of this Figure again reveals that the maximum disturbance point 
does, in this case, exist at a tangential location that is lagging the minimum-clearance point 
(relative to the impeller speed). It should be emphasized, however, that our objective here 
was no more than a qualitative assessment of our results, particularly in view of the fact 
that no pressure measurements were reported for a whirling seal rotor neither by Morrison 
nor by any other experimentalist to the author’s best knowledge. 


Validation of the Single-Harmonic Perturbation Assumption: 

The assumption that perturbations in each flow property is expressible via a single- 
harmonic tangential distribution, has commonly been used in all perturbation models 
existing today [e.g. Childs (1989), Dietzen and Nordmann (1987),...]. The current pertur- 
bation model, on the contrary, does not pre-impose any such restriction. 


In order to investigate the validity of this assumption, the circumferential distribution 
of the shroud pressure perturbation was expressed at each through-flow location (using 
the already-computed nodal values) as a Fourier series as follows: 


dp ^ 

~ a ° + ^ cos (nO) + ^ b n sin(n0) 


where “0” is the angular coordinate. The number of sine and cosine waves “n” in this 
expression was arbitrarily set equal to three. The ratios ( a n ja \ ) and ( b n /b \ ) were then 
examined for n = 2&3. The outcome of this study was that these ratios were sufficiently 
close to zero (meaning a valid single-harmonic assumption) at locations away from the 
leakage-passage inlet and exit stations and away from the passage segments where excessive 
flow recirculation and vortex breakdown occur (Fig. 8). In fact, the ratios a n j a\ and b n /b\, 
where n = 2&3, can be as high, in the latter subdomains, as 21%. Fig. 17 shows the shroud 
segments where the single-harmonic assumption is hardly accurate. 








PROPOSED WORK 


Falling under the category of rotordynamics, three different expansion studies are 
proposed. Of these, the first proposed study concerns creation of a new and desirable 
mode of code execution, whereby the phase of solving the centered-rotor flow field is 
by-passed, and a pre-existing flow solution used instead. The second study concerns an 
upgrade of the perturbation model to account for contributions of the primary flow to the 
shroud forces. The third proposed study involves what promises to be the most rigorous 
rotordynamic analysis of hydrostatic bearings. A summary of the objectives and approach 
in each study is given hereafter. 


1. Insertion of an Existing Unperturbed Flow Field: 

As described in Part 2 of this report, a new option (in executing the computational 
model) is proposed, whereby the user may make use of an existing centered-rotor (ax- 
isymmetric) flow solution [or, more importantly, centered-rotor flow measurements (if 
available)] as input to the perturbation model. The user, in this case, would be overriding 
the subprogram which is responsible for securing this “zeroth-order” flow solution. Due 
to the vast differences in the grid structure and the computed (or measured) flow variables 
there can be (let alone the format of the user’s files), the author realizes the difficulty 
in “universalizing” this option. The focus will therefore be on providing the user with 
a systematic procedure whereby he/she could write his/her own interface code (usually 
short and simple), which would convert the grid (and, occasionally, the flow variables) into 
the form and format that is acceptable to the current perturbation analysis subprogram. 
In addition, at least one example of such interface code, the present outcome of which 
is the middle grid in Fig. 18, will be provided in the final version of the User’s Manual. 
Nevertheless, it is important to point out that the user would (from a practical stand- 
point) resort to this option only if he/she has a good reason to (e.g. availability of reliable 
experimental data, or a highly accurate axisymmetric flow solver featuring, for instance, 
a state-of-the-art turbulence closure, ...etc.). Otherwise, the user is generally advised to 
use the existing centered-rotor flow solver. 

Included in this first proposed study is also the execution of a bench-mark test case, 
which was recently reported by Cal. Tech. (Guinzburg, 1992). The significance of Cal. 
Tech.’s rotordynamic measurements stems from the fact that the rig was designed to ex- 
clusively account for the leakage flow effects only, meaning that the primary (impeller) flow 
in this rig was non-existing (Fig. 19). As reported by Cal. Tech., the fluid-exerted forces 
on the rotor were parabolic and fluctuation-free, when plotted versus the whirl frequency. 
The author believes that this test case will complement the ongoing code verification pro- 
cess, and have actually started the rig simulation process by creating the centered-rotor 
finite element grid (Fig. 19). 
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Inclusion of the Primary Flow Passage in the Analysis: 

The current version of the currentperturbation model only accounts for the shroud 
forces developed in the secondary (leakage) passage (Fig. la). Since the inner shroud 
surface is also exposed to the flowing medium, it would be inaccurate to extract the latter 
from the flow domain with which the shroud interacts. 

We are proposing a two-phase study to account for the primary flow passage effects 
on the shrouded-impeller rotordynamic characteristics. Referring to Fig. 20, these two 
phases are as follows: 

a) Full guidedness of the impeller flow will first be assumed, whereby the swirl velocity at 
any finite-element node within the primary-flow passage (Fig. 20) will be determined from 
an already computed meridional (r-z) velocity and the assumption that the total velocity 
vector is tangent to the impeller blade. This is an iterative procedure the outcome of w T hich 
will be combined with the leakage-passage flow field prior to executing the perturbation 
model. 

b) A new computational capability, which we have just developed, will be utilized for a 
more realistic simulation of the primary-passage flow field. This capability is a blade- 
to-blade turbulent flow 7 analysis code for a centrifugal pump impeller. The code is now 
operational, and a sample pump impeller case for w'hich experimental data, by Lanes 
(1982), was provided. This study has already been concluded (Hlavaty, 1993) under the 
author’s supervision. Results of this code, for a given impeller geometry and operating 
conditions, will be used to compute the circumferentially-averaged swirl velocity compo- 
nent which, in turn, will will be used in the axisymmetric-flow governing equations to 
compute a comparatively accurate zeroth-order flow solution in the primary— flow passage. 
This solution will then be merged with the leakage-passage axisymmetric flow solution, 
with the outcome being a valid input to the perturbation model. 

As indicated above, we do have the computational means of implementing these lower- 
and higher-order modeling phases. In the end, we will be in a position to assess the 
contribution of the combined primary/secondary flow passages, in a real pump situation, 
versus the contribution of only the leakage-flow passage. 


Rotordynamic Characteristics of Hydrostatic Bearings: 

This is perhaps the most aggressive, yet very much achievable, proposed task. The 
proposed approach (summarized hereafter) is currently under consideration, and should 
cover a range of bearing between outlined hereafter is applicable to a range of bearings be- 
tween totally pressurized (hydrostatic) to hydrodynamic bearings including, in particular, 
the hybrid bearing category. In our model, there will be no assumption of viscosity- 
dominated (low Reynolds’ number) flow anywhere, since the inertia terms will always be 
retained, no compromise of the flow turbulence anywhere (and not only over solid walls), 
and no compromise of the journal-to-housing gradients of the flow variables. In the sense 
that 
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Fig. 20 Inclusion of the primary flow passage as part of the computational domain 




the model is a natural extension of what we generically termed the “virtually” deformable 
finite-element approach, it should also be emphasized that utilization of the finite-element 
method in the current analysis is not only meant to solve one or more orders of the flow 
perturbation equations, but is also the means with which these equations come to existence. 
This is true in the sense that the perturbation equations, in the current model, emerge 
from expansion of the flow-governing equations in their discrete finite-element form, as 
opposed to the differential form in any existing perturbation model. This feature is the 
reason behind the model unmatched versatility, especially when it comes to the three- 
dimensionality of the flow domain. 

In the proposed model, a time-accurate algorithm is utilized to determine the history 
of the lubricant-exerted forces on the journal, as the latter traverses over a tangential 
“pitch” between two successive recesses. This will make it possible to compute a set of 
rotordynamic coefficients, which are based on the time-averaged fluid-exerted forces, and 
another set of “conservative” coefficients which correspond to the instance at which the 
destabilizing/stabilizing journal-force ratio is maximum. In the following, the proposed 
approach will be summarized with reference to Fig. 21, which shows the NLS Fuel Pump 
hydrostatic bearing as an example. 

First, the reader is reminded of some relevant details concerning the current pertur- 
bation model. Referring to Fig. 22, the frame of reference in which the “unperturbed” 
flow field is described consists of the z-axis (which at this point is coincident with the 
housing centerline) with the x and y axes rotating at an angular speed that is identical to 
the whirl frequency fl. Recalling that the journal whirl has not physically occurred yet, it 
is seen that this coordinate system is highly untraditional but is, nevertheless, legitimate. 
The reason for this choice is that once the governing equations and boundary conditions 
are cast in this coordinate system, the only physical component that is missing from the 
whirling motion is the journal eccentricity e. In other words, if the journal whirl is thought 
of as a journal eccentricity coupled with a whirl frequency, then the latter is taken into 
account in the early step of defining the centered-journal flow field. Furthermore, because 
the flow variables (including the velocity) are relative to this coordinate system, it is actu- 
ally the housing which undergoes the rotational motion at this point. This is simply due 
to the fact that the flow variables and boundary conditions are those which would appear 
to an observer in this rotating frame of reference. Fig. 6 shows the solid-wall boundary 
conditions relative to that observer. 

With the preceding logic in mind, and referring to the hydrostatic bearing problem (Fig. 
22), all one has to do for determining the perturbations in the fluid properties (particularly 
those of the journal-surface pressure) is to construct a centered-journal (zeroth-order) flow 
problem that is consistent with that associated with the position of the whirling journal 
at any given point in time. The consistency here includes the instantaneous orientation 
of the coordinate axes and the whirl frequency Q as shown in the upper row of Fig. 22. 
Superimposition of the journal lateral eccentricity e (lower row of Fig. 22) would then 
complete the requirements for the whirling motion (at the corresponding point in time) to 
prevail. The procedure is therefore one of marching in time, in which the journal 
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Fig. 22 Instantaneous positions of the coordinate axes over one 
angular pitch between two successive recesses 
in a hydrostatic bearing 



repeatedly sweeps (in its whirling motion) the angular spacing between two successive 
recesses (Fig. 22), until fixed cyclic patterns of the fluid reaction forces are achieved. 
During the marching process, a “centered-journal” flow field is first obtained for any given 
angular position of the journal axis (upper row in Fig. 22). Subsequent application of 
the “virtually” deformable finite element concept (Baskharone and Hensel, 1991a) to this 
zeroth-order flow field will then produce the perturbations of the flow variables. Of these, 
the journal-surface pressure perturbations are extracted and integrated over the journal 
surface to attain the instantaneous radial and tangential components (F r and Fq ) of the 
journal/fluid interaction force. 

The procedure described above will produce the history of the fluid-reaction force com- 
ponents as the journal whirls through a complete angular “pitch” between two successive 
recesses, which is 60 deg. in Fig. 22. The time-averaged values of these force components 
can then be computed and used in a traditional rotordynamic analysis, as outlined earlier. 
Note that the proposed approach compromises none of the real flow effects, and fully ac- 
counts for the three-dimensionality of the journal-to-housing domain. It is also important 
to point out that despite the generality and rigorousness of the proposed approach, the size 
of the model (in terms of computational resources) is not expected to be much different 
from that of the shrouded-impeller rotordynamics problem, which we have successfully 
solved . 
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A Finite-Element Perturbation 
Approach to Fluid/Rotor 
Interaction in Turbomachinery 
Elements. Part 1: Theory 

The vibrational characteristics of a rotor that is in contact with a fluid in an annular 
clearance gap , as dictated by the fluid forces in the gap, are investigated. The “rotor” 
here is a general term that may refer to the shaft segment within the housing of an 
annular seal \ on the simple end of the application spectrum, or the shroud-seal 
assembly in a shrouded-impeller stage of a turbomachine, on the complex end. The 
disturbance under consideration involves the axis of rotation, and includes a virtual 
lateral eccentricity, together with a whirling motion around the housing centerline . 
Uniqueness of the computational model stems from the manner in which the rotor 
eccentricity is physically perceived and subsequently incorporated. It is first estab- 
lished that the fluid reaction components arise from infinitesimally small defor- 
mations with varied magnitudes which are experienced by an assembly of finite 
elements in the rotor-to-housing gap as the gap becomes distorted due to the rotor 
virtual eccentricity . The idea is then cast into a perturbation model in which the 
perturbation equations emerge from the flow-governing equations in their discrete 
finite-element form as opposed to the differential form, which is traditionally the 
case. As a result , restrictions on the rotor-to-housing gap geometry, or the manner 
in which the rotor virtual eccentricity occurs are practically nonexisting. While the 
emphasis in this paper is on the theoretical model, a representative application of 
the model and assessment of the numerical results are the focus of a companion 
paper that is being published concurrently. 


Introduction 

Reliable modeling of the rotordynamic characteristics of 
fluid-encompassed elements in turbomachines, such as the an- 
nular seal rotor in Fig. 1, started fairly recently. The bulk- 
flow model devised by Childs (1983) is an example of some 
widely used predictive tools in this area. This technique was 
used to compute both the direct and cross-coupled rotordy- 
namic coefficients for an annular seal (Childs, 1985) and fur- 
ther extended by Nelson (1985) to account for the flow com- 
pressibility and the seal taper. However, the accuracy of this 
model would naturally decline in the case of a generally shaped 
annular flow passage, such as the secondary passage of the 
shrouded impeller in Fig. 2, where lateral (rotor- to-housing) 
flow stresses become comparable to those at the solid walls. 
Besides, the analysis is conceptually incapable of simulating 
the effects of flow separation and recirculation in such a com- 
plex flow passage, but has been remarkably successful oth- 
erwise (Childs, 1989). 


’Graduate Research Assistant, currently a Senior Engineer at Westinghouse 
Savannah River Laboratory, Aiken, South Carolina. 
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Recently, a more rigorous Finite difference-based pertur- 
bation model was developed by Dietzen et al. (1987) and ap- 
plied to an annular seal. The distorted seal annulus for the 
displaced shaft configuration, in this case, was mapped into 
a Fictitious frame of reference, reducing the problem into that 
of a centered shaft rotation. A traditional hypothesis was then 
made that perturbations of the flow variables vary sinusoidally 
with the tangential coordinate around the circumference in the 
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Fig. 2 Two representative problems of turbomachinery rotating elements 
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fictitious space, with only the first harmonic wave taken into 
account. Based on this assumption, a perturbation analysis of 
Navier-Stokes equations was then performed and the rotor- 
dynamic coefficients, in the end, obtained. Despite the ap- 
parent generality and complex nature of this analysis, it involves 
geometry-related restrictions and is limited to a uniform lateral 
eccentricity of the rotor axis. Other existing models for annular 
and labyrinth seals, e.g., those by Dietzen and Nordmann 
(1988) and Hensel (1986), do not utilize the perturbation con- 
cept and are therefore tedious as well as costly, for they require 
solution of the fully three-dimensional Navier-Stokes equa- 
tions in the distorted flow domain. 

In this paper, a new finite element-based perturbation model 
is presented. The model introduces what may be referred to 
as the “virtually deformable” finite element concept, as a 
means of determining the fluid-exerted forces in the rotor-to- 
housing gap, and whether such forces would sustain or suppress 
a rotor whirling motion of the cylindrical type, as depicted in 
Fig. 1 for a simple annular seal. Originality and versatility of 
the new model both stem from the fact that the perturbation 
equations emerge from the discrete finite-element flow equa- 
tions as opposed to the flow-governing equations in their dif- 
ferential form as is traditionally the case. 

Analysis 

The governing equations are grouped as pertinent to the 
fundamental rotordynamics basis, the flow kinetics in the un- 
disturbed operation mode and, finally, the perturbation model. 
In all phases, reference is made to a general fluid domain with 


which the rotor is in contact. Two examples of such domain, 
namely that of an annular seal and shrouded pump impeller 
(Fig. 2), are given. 

Definition of tbe Rotordynamic Coefficients. Evaluation 
of the force coefficients is initialized by considering the rotor 
lateral eccentricity “c M coupled with a circular whirling motion 
around the housing centerline (Fig. 1). The locus of the rotor 
axis, in this case, can parametrically be described as follows: 

X ~ esin(ftf ) , Y = ecos(G/ ) ( 1 ) 

where “X” and “7” are the eccentricity components of the 
rotor axis (Fig. 1), “G” is the whirl frequency and “r” refers 
to time. Noting that the direction of “0” in Fig. 3(a) is neg- 
ative, the rotor equation of motion can be written as follows: 



where 

bF xy bF y are the incremental changes in the rotor force com- 
ponents as a result of the rotor disturbance, 

AT, k are the direct and cross-coupled stiffness coefficients 
of the fluid/structure system, 

C, c are the damping coefficients, 

A/, m are the inertia (or added mass) coefficients, 

Assuming cylindrical whirl around the housing centerline at 




V 


Nomenclature 


[A], [a] 

IB), ib) 
F, G, H, I, J, K 
K y k t C, c, Af, m 


p 

Uy Vy W 


W i 


global and elemental matrices of in- 
fluence coefficients in the finite ele- 
ment equations 

global and elemental load vectors in 
the finite element equations 
linear operators defined in the element 
local frame of reference 
direct and cross-coupled stiffness, 
damping and inertia coefficients 
linear shape function associated with 
the rth corner node of a finite element 
quadratic shape function associated 
with the rth corner or midside node of 
a finite element 
static pressure 

velocity components in the cartesian 
frame of reference 
weight function in the Petrov-Galer- 
kin weighted-residual analysis 


y, z 
v> £ 

v 

1 * 1 . 1*1 
Q 

co 


Other Symbols 

overbar 

O 


cartesian coordinates 

lateral eccentricity of the rotor axis 

coordinates in the element local frame 

of reference 

kinematic viscosity 

density 

global and elemental vector of un- 
known velocity components and pres- 
sure 

whirl frequency 
rotor operating speed 


perturbation in a quantity due to the 
rotor eccentricity 

a quantity that is known from a pre- 
vious iteration or an initial guess 
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Fig. 3 Definition of the rotating-translating axes and solid-wall 
boundary conditions 


a constant whirl frequency M fl’\ and refering to Fig. 1 , consider 
the position of the rotor axis that is vertically underneath the 
housing centerline. The axis location, linear velocity and ac- 
V celeration at this position are: 

X=0 Y=t 
X=Qe Y= 0 

x=o r=-n 2 e 

Substituting these into the expanded form of the matrix Eq. 
(2), and taking the limit as ‘V* tends to zero, the following 



equations are obtained: 




bF x /& 

F x \ , 




=k-nc-tfm 

(3) 


dF, (6F 

-r z = lim( 
be ( -o \ e 

A = -K-Qc + n 2 M 

(4) 


V Determination of the rotordynamic coefficients in Eqs. (3) and 

dF bF 

(4) requires computation of the derivatives — - and — * at a 

be be 

minimum of three different values of the whirl frequency “0.” 
Interpolation of these two derivatives as parabolic expressions 
of “fl,” using curve fitting techniques, leads to the rotor- 
dynamic coefficients by simply equating the different terms in 
these expressions to those on the right-hand sides of Eqs. (3) 
and (4). It is therefore clear that the fundamental problem 
here is that of computing the rate (with respect to the virtual 
eccentricity “e”) at which the fluid reaction forces are devel- 
oped on the whirling rotor (Fig. 1). The procedure is initialized 
by solving the flow-governing equations in the undisturbed 
rotation mode as presented next. 

v Choice of the Coordinate Axes. The flow equations are 

cast in a frame of reference that is consistent with the annular 
gap geometry in the disturbed-rotor operation mode. As shown 
in Fig. 3(a), the reference axes are attached to the displaced 
rotor and are, therefore, of the rotating-translating, or simply 
whirling, type. The rotational speed, in this case, is numerically 
equal to the whirl frequency “O'* and the translation is identical 

V to that of the rotor axis. As seen, the clearance gap flow field, 
viewed in this frame of reference, is now steady at all times, 
despite the apparent time dependency of the absolute flow 
field when viewed in a stationary frame of reference. 

Motion of the coordinate axes, in the general manner just 
described, takes place only after the rotor has entered its ec- 
centric operation mode. However, it is only the axes rotation 
that prevails in the centered-rotor mode (Fig. 3(6)), in which 
case “0" is to be viewed as an arbitrary frequency. 

The choice of a rotating-translating frame of reference does 
not imply the need to solve the flow-governing equations in 
the physically distorted clearance gap (e.g., Dietzen et al., 
(1987). The fact, however, is that this choice, despite the com- 


plex nature of the flow equations it produces (as will be seen 
next), makes it easier to develop the perturbation model as it 
ensures consistency of the unperturbed and perturbed flow 
fields, with the difference between the two Fields being a result 
of only the rotor eccentricity (Fig. 3). Note that this consistency 
requirement can alternately be ensured by solving the simple 
axisymmetric flow field in the undistorted passage (Baskharone 
and Hensel, 1991) and then casting this “zeroth-order” so- 
lution in the frame of reference defined above (Hensel, 1990). 
This procedure was implemented during execution of the com- 
puter code for the substantial reduction in the core and CPU 
time consumption it provides. 


Flow-Governing Equations. The swirling flow in the un- 
distorted rotor-to-housing gap (Fig. 3(6)) is assumed adiabatic, 
incompressible and generally turbulent. The momentum and 
mass conservation equations can therefore be expressed in the 
rotating frame of reference (noting that “0” in Fig. 3(6) is 
negative) as follows: 


. du n du „ du ^ 1 bp 2 

tf T -+ 0 — + w — + 2Ql> - Q 2 x= - - — + v e V l u 
bx by dz p bx 


bv ( du bv t (bu 3iA bv t (bu dw\ 
bx bx + by \dy + bx) + bz \dz + bx) 


„ bv „ bv „ bv ^ 1 bp 2 

u T + v 7“+ w - - — + v e V v 

bx by bz p bx 


bv t bv biff/bu du\ bv ( /bv 
+ by by* bx \dy + bx) + bz \dz + by) 


„ bw „bw „ bw 1 bp 2 

a —+o — + w — +^v 2 w 

bx by bz p bz 


+ 2 


bv t bw bv t (bu dw\ 
bz bz + bx + bx) 


bu bv bw . 

T" + T~ + — - 0 

bx by bz 


(5) 


( 6 ) 


(7) 

( 8 ) 


where 

w, u, w are the relative velocity components 
p is the static pressure 

0 is the rotational frequency of the coordinate axes (to 
be set equal to the rotor whirl frequency once the rotor 
eccentricity prevails) 

v tt v e are the eddy and effective coefficients of kinematic 
viscosity, respectively 

with the symbol O designating values that are carried over 
from a previous iteration or an initial guess for the purpose 
of successively linearizing the momentum equations. As seen, 
the axes rotational frequency “Q” is part of these equations, 
which is a result of including the centripetal and Coriolis ac- 
celeration effects in this case. 
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Turbulence Closure. Simulation of the flow turbulence is 
based on the algebraic eddy-viscosity turbulence model by 
Baldwin and Lomax (1978). However, an adjustment in this 
vorticity-based model was implemented whereby the relative, 
as opposed to the absolute, vorticity was used to calculate the 
mixing length in the fluid layer that is adjacent to a solid 
surface. Equally important in implementing the turbulence 
model was the analysis of the near-wall zone which was con- 
ceptually based on the approach by Benim and Zinser (1985). 
Further details of the turbulence closure and the near-wall flow 
analysis are covered by Baskharone and Hensel (1991). 

Boundary Conditions. The major difference between the 
two problem configurations in Fig. 2 is in the multiple entry/ 



a) LOCAL FRAME OF REFERENCE 



Fig. 4 Quadrilateral curve-sided isoparametric finite element 


departure nature of the flow domain associated with the cen- 
trifugal impeller case. Uniqueness of the solution, under such 
circumstances, was described and attained by Baskharone and 
Hensel (1989) through a specific set of boundary conditions 
over the flow-permeable segments of the domain boundary. 
Boundary conditions pertaining to the annular seal problem 
include known profiles of inlet velocity components and zero 
streamwise diffusion of the exit velocity vector. The latter 
replaces a zero velocity gradient exit condition which would 
imply a fully developed flow at this station and may, in this 
sense, be accurate only for sufficiently long seals. 

The last category of boundary conditions concerns the rotor 
and housing surfaces and is created by the rotation of the 
coordinate axes. Referring to Fig. 3(£), the housing surface 
to an observer in the rotating frame of reference, will no longer 
appear stationary, but will rather possess a relative velocity 
component “fib” as indicated in the figure, where il b ” is the 
housing radius. Considering the case of what would amount 
to forward whirl (as the rotor axis becomes eccentric), the 
rotor surface in Fig. 3(6) would appear to the same observer 
as rotating at an angular speed that is less than the rotor 
operating speed by the amount “0.” The rotor surface 
velocity in this case is (u> - ffyz, where "a” is the rotor radius. 

Finite-Element Formulation. The finite-element model is 
constructed with the twenty-node quadrilateral element in Fig. 
4 as the discretization unit. The procedure is initiated by cre- 
ating the finite-element discretization model corresponding to 
the centered-rotor operation mode as shown in Fig. 5(a) for 
a generally-shaped rotor. Within a typical “undistorted” finite 
element in this figure, let the cartesian-to-local spatial coor- 
dinate mapping (Fig. 4) be defined as follows: 

20 20 20 

x=J]N i x i , ( 9 ) 

i»l /«1 im\ 

where N,- are quadratic “shape” functions (Zienkiewicz, 1971) 
of the local coordinates f, tj, and £ (Fig. 4). Conversion of the 
spatial derivatives is, in this case, defined as follows: 




Fig. 5 Distortion of the finite elements cs a result of the rotor eccen- 
tricity in s general fluid/rotor interaction problem 
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where the components of the matrix [7] are as follows: 



with 171 being the Jacobian of transformation. 

Next the flow variables are interpolated within the typical 
element in a similar fashion. Guided by the Ladyshenskaya- 
Babuska-Brezzi compatibility requirements (Carey and Oden, 
1986), as applied to the current problem, the velocity com- 
ponents and pressure are expressed as follows: 

20 20 

u = 2 N ‘ u « v = 2 N ' v ” 

r«l /-I 

20 8 

w= NjW h p = 2 M*Pk (11) 

fT\ *-i 

where the interpolants “M*” appearing in the pressure expres- 
sions are the linear shape functions associated with the finite 
element corner nodes (Zienkiewwicz, 1971). 

The error functions produced by the unperturbed flow gov- 
erning equations as a result of the interpolation expressions 
above, are then made orthogonal to a special set of weight 

functions throughout the finite element. In constructing 

these functions, the so-called error consistency criterion of 
Hood and Taylor (1974) is implemented, whereby the element 
linear shape functions, Af;, are used in conjunction with the 
continuity equation. The weight functions used with the mo- 
mentum equations, on the other hand, are of the upwind form 
(Baskharone and Hensel, 1991) when applied to the convection 
terms, and are identical to the element quadratic shape func- 
tions for all other terms. 

The orthogonality conditions stated above constitute the 
elemental set of finite element equations, which can be written 
as follows: 



+ l«dlwl + folM + l««llPl = f*il 

(12) 


(a 5 ]{u) + [a 6 ] 1 v ) +[ 07 ](h’] +[a 8 ]|p) = (* 2 ) 

(13) 


[fl,](w) +[a| 0 ]{t'i + [«n](w) + [fli2](/>l = (* 3 ) 

(H) 


[fli3lf«) + [0u](u) +(<ri5l(w) =0 

(15) 

where: 




( (P,[FW)F(^) + G(Ni)G(Nj) + H(N,)H(N > )] 

Jk(») 


-Nj[G(v t )G(Nj) + H(»,)H (Nj) + 2F(P / )F(A r y )] 

4- Wj[uY(Nj) + 0G(Nj) + wHINj)) )dV 

a 2 ) itJ = [ -*,[G(P,)F(A0) -2UNjldV 

J y(e) 

[ -NmnNj)dV 
Jyie) 

a 4 ) ; .*=( - N,F (M k )dV 

Jk(') P 

a s ),j= ( - M[F(P,)G (Nj) + 2U Nj]dV 

<dy[t) 

<J 6 ), j=\ ( ?,[F(N,)F (Nj) + G (Ni)G(Nj) + H(/V,)H(N,)] 

-yV,[F(P r )F(N,) + H(*,)H(ty) + 2G(P,)G(N,)] 

+ ^,[i3F(N y ) + vG(Nj) + wH (Nj)] \dV 

),.,= ( - NiH(v,)G(Nj)dV 

~NG(M k )dV 

Jy[e) P 

«»)/.;= ( - N,¥(v,)H (Nj)dV 

«,o),.;= ( -A/,G(i>,)H(N,)</K 

Jk(») 

a„),.,= f [v t \¥(.Nj)T(Nj) +G(Nj)G(Nj) + H(N,)H(Ny)] 


-N,[F(P,)F(N;) + G(v,)G(Nj) + 2H(P,)H(N,)] 

+ W,[uY{Nj) + vG(Nj) + wH(N ,)] \dV 

-N,H(M t )dV 

J y(t) P 

a i3 ) kJ = ( (Nj)dV 

J y(e) 

a„) kJ = \ M k G(Nj)dV 

* 5 )*.,= ( M k H(Nj)dV 

J y(t) 

bi),= ( 

Js<«) \/»i / 

62),= ( v t N i (n-W)dS + O 2 f N,(2 N '> , 'r K 

Js<«> j k<*> \i - 1 / 

f> 3 )/= \ v l N l (n-'7w)dS 

where ( e ) refers to a typical finite element, / = 1,2 ,..., 20, 

jm 1,2 20, k = 1, 2 , . . . , 8, the superscript O refers 

to a value that is known from a previous iteration or an initial 
guess, dS is the differential element of surface area and n is 
the local unit vector normal to the boundary (Fig. 4), Also, 
the volume differential dV and the operators F, G, and H are 
defined as follows: 

dV= \J\d?dvd£ 


V 
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H - r »r r -4 +7 ' ! 4 ,l8) 

It should be pointed out that expanding the finite element 
equations in terms of the operators F, G, and H, defined above, 
is consistent with, and indeed simplifies, the procedure leading 
to the perturbed version of the flow equations, as will be seen 
later in this section. 

Equations (12) through (15) can be rewritten in the following 
compact form: 

= (19) 

where the vector (<£) contains the nodal values of velocity 
components and static pressure that are associated with the 
typical element “e.” The global form of the elemental set of 
Eqs. (19) is obtained by assembling all finite element contri- 
butions, and introducing the various boundary conditions, at 
which point the equations are expressible in the following form: 

H](*) = {*) (20) 

Perturbation Model. Referring to Eqs. (3) and (4), deter- 
mination of the rotordynamic coefficients reduces to the fun- 
damental problem of computing the rate (with respect to the 
virtual eccentricity 4 V) at which the fluid forces are exerted 
on a cylindrically-whirling rotor. This objective is achieved 
through the perturbation analysis outlined next. 

The current model is centered around the manner in which 
the Finite-element Eqs. (20) are altered as a result of the rotor 
virtual eccentricity “e” (Fig. 5). Under such eccentricity, each 
finite element will yield a perturbed set of equations. The 
assembled form of these equations can generally be written as 
follows: 

(lA] + elA])(l*} + l6*)) = (lB)+tlB)) (21) 

where ($) is the global vector containing the nodal values of 
the velocity components and pressure. Combination of Eqs. 
(20) and (21) gives rise to the following: 

3($) ( 6 *) , - - 

J=M]- , aB}-H]|*)) (22) 

Equation (22) reveals that the differential changes by which 
the flow variables, namely the velocity components and pres- 
sure, vary as a result of the rotor eccentricity can be achieved 
knowing the undisturbed flow solution vector ( $ ) which is 
known at this point, a matrix [A] and a vector (2?) with the 
latter two arrays representing, respectively, the effect of dis- 
torting the finite elements, and the changes in the flow kine- 
matics as a result of the whirling motion of the coordinate 
axes (Fig. 3(a)). The arrays [A] and (£) are derived next as 
part of the procedure to compute the rotordynamic coeffi- 
cients. 

Consider a distorted finite element in the clearance gap cor- 
responding to the eccentric-rotor operation mode (Fig. 5). To 
an observer in the whirling frame of reference, the typical node 
of this element will now be displaced by an amount that 
is a function of the rotor virtual eccentricity ‘V* and the node 
original location in the cross-flow plane. Referred to the dis- 
placed axes position in Fig. 5(2?), the new nodal coordinates 
x„ y„ and z , can be related to those prior to the rotor eccentricity 
as follows: 


X/ — Xj f y> — y i + X/€ and Z; — z ( 

where X, is a fraction that varies from zero, for nodes on the 
rotor surface, to — 1.0 for those on the housing surface, since 
the housing, to an observer in the whirling frame of reference, 
is the surface that undergoes the virtual displacement (note the 
positive direction of the .y-axis in Fig. 4). The local-to-cartesian 


transformation Jacobian I 71 can correspondingly be written 
for a distorted finite element as 1 7, 1 where: 

17,1 = 17 1 +c 171 (24) 


in which the Jacobian I 71 is that of the undistorted element, 
and I 71 is as follows: 


171 = 


20 dN, ^ dN, 


20 


V 2fZi v V °' V| •> v> dN < 


20 


V x V — X V ^ 


20 dN, dN, 


E ojy, V \ V 1 dN ‘ 

. ^ di Xi dk z ' 


A new spatial derivative operator, that is equivalent to that in 
Eq. (10), can now be derived as follows: 



where the matrix [7] is the same as described in conjunction 
with Eq. (10), while the matrix [/>] is defined as follows: 
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Using the Jacobian of coordinate transformation and operator 
of spatial derivatives for the distorted element (Eqs. (24) and 

(25) , respectively) as well as the general interpolation expres- 
sions (1 1), in the flow-governing Eqs. (5 through 8), and reap- 
plying the weighted-residual procedure, the distorted element 
equations can be written (upon separation of terms containing 
4 V* and ignoring high order terms) as follows: 

(w+€ra)([*i + i««!)®(i*i +€[*)) (26) 

where the matrix (5) and the vector ( b ) are produced by the 
distortion in thefinite element shape, while the non-zero entries 
of the vector |Z?) are associated with some distortion-related 
entries which are explained later in this section. The vector 
( b) also contains the effect of an additional acceleration com- 
ponent that is created by the motion of the coordinate axes 
origin around the housing centerline. The vector { 6<t> } in Eq. 

(26) contains the differential changes in the nodal values of 
velocity components and pressure that are created by the virtual 
eccentricity. 

Noting that the matrix [a] is similar, in construction, to [a] 
in Eq. (19), the earlier can consistently be viewed as composed 
of submatrices [ari through which correspond to those 
in Eqs. (12) through (15), and are defined as follows: 

= *?),,,+ ( (^IF(M)I(N,) 

J y(t) 

+ m)F(N J )+G(N i )Wj) 

+ i(N,)G(N J )+U(N i )K(Nj)+K(N l )H(Nj)]-N l \G(t',)3(N J ) 
+ J(P,)G (Nj) + H(P,)K(N,) + K(P,)H(Ny) 

+ 2F(P,)I(AT,) + 2I(P,)F(N,)] + WmNj) 

+ \i(N J ) + w¥.(N J )])dV 

<i2hj = alhj- [ M[G(P,)I(N,)+J(P,)F(Ny)]tfK 

J y(e) 

a 3 )u=«*)u- [ N<[H(P r )I(Afy) + K(P,)F (Nj)]dV 

Vy(t) 

5«),.* = a;) ( ,*+ [ -N,l(M k )dV 

Jk(»> P 

a s )ij = a? ) u - ( Af,[F(P,)J (Nj) + I(P,)G (Nj)]dV 
J(/«) 

Ot)i,j = a i)ij+ ( (P,[F(/Vi)I(N,) 

Jk(») 

+ l(N)F(N J ) + G(N i )3(N J ) 

+J(N,)G(Nj)+H(N i )K(Nj)+K(N i )H(N j )]-Nmh)HNj) 
+ I(P,)F (Nj) + H(P,)K(N,) + K(P,)H(A r y ) 

+ 2G(P,)J (Nj) + 2J(P,)G (Nj)] + WmNj) 

+ 03(Nj) + wK(Nj)]]dV 

d-,) u = a f ( Ni\H(v,)J(Nj) + KC»,)G(Nj)]dV 

J y(e) 

= ( -N^mdV 

J y{e) P 

= [ Ni[Fft)K(N,)+Ift)H(Af y )]£/K 

Jyie) 

5,o) (V = afo)i.>- ( M[G(P,)K(N,)+J(P,)H(N,)JdK 

J y{e ) 

«.i>w-«?.)u+ ( (P,[F(N,)I(N y ) 

+ \(Ni)¥ (Nj) +G(Nj)J(Nj) 


+ J(N,)G(N i ) +U(N t )K(Nj) ^KiN^HiN^-NmW^) 
+ 1( h)F(Nj) + G(y,)J (Nj) +J(v t )G(Nj) 

+ 2H(J r )K(/^) +2K(P,)H(Af y )] + W t (u\(Nj) 

+ vJ(Nj) + wK(Nj)]\dV 


a \2)i,k = <*\2)i,k + \ 

- N,K(M k )dV 

J 

P 

Q\l)kJ = <*h)kJ + 

\ M k \(Nj)dV 


•V<»> 

Q\4)kJ = a U)k,j + 

( M k 3(N J )dV 


Jk(') 

a\5)kj = < 1 l$)kj + 

[ M k K(Nj)dV 

Jy<e) 


where 

dV= \J\dtfridZ 

In the these expressions, a matrix entry with an asterisk has 
the same form as that in Eqs. (12) through (15), with the 
exception that the Jacobian 171 is now replaced by 171 as 
defined in conjunction with Eq. (24). The vector |£) in_Eq. 
(26) can similarly be defined by the subvectors ( b\ ) and ( b 2 ] , 
as in Eqs. (12) through (15), as follows: 


S,),-*^ Nfi \J\d{di,dk 

' = 0 l L.l N, (s N °^) 

-tff N/dV 

J ^'> \Jt| / «V<»> 


*3)/ = 0 

where i and j vary from 1 to 20, while k varies from 1 to 8. 
With the operator F, G, and H being those defined by Eqs. 
(16) through (18), the new operators I, J, and K are defined 
as follows: 


d d d 


I=p,1 ar +p,2 d^ +p ' J a£ 

(27) 

d d d 

J = p 2l + 

(28) 

a d d 

K=p "K +p »T v +Pli K 

(29) 


The matrix [P] is the same as defined in conjunction with Eq. 
(25)._Of particular interest here is the construction of b\) n b 2 )( 
and by)i above, as compared to that of b])*, b 2 )j and b 3 )/ in 
Eqs. (12) through (14). First, it is noted that the perturbations 
of the surface integrals in bj), and Z> 2 ), are nonexisting in b { 
and b 2 since these surface integrals would have nonzero values 
only in the case where a nonzero normal derivative of, re- 
spectively, ‘V’ or * V* is prescribed as a boundary condition, 
a situation which is not encountered in formulating the un- 
disturbed-flow problem. Secondly, it is due to the centripetal 
acceleration of the rotating axes origin around the housing 
centerline that the last term in the b 2 )> is produced, while the 
other terms, in b\)i and b 2 ) h result from infinitesimally-small 
distortions in the finite element shape as a result of the rotor 
eccentricity. 

Calculation of the Rotordynamic Coefficients. Equation 
(22) can be rewritten in the following detailed form: 
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be 


d_ 

be 


( {U] ) 

m 

\m 

U P)J 


= [A]- l UB}-[A][*}) 


(30) 


where ((/), { K) , [W) t and \P] are vectors containing the 
nodal values of the velocity components and static pressure, 
respectively, with the overbars signifying quantities that are 
associated with the perturbation in the rotor-to-housing flow 
Held. Of all the pressure nodal values in the vector (P) above, 
define a subvector ((P) as composed of the pressure values at 
the rotor surface nodes, i.e., 


\<P} = lp,J=l,N s }C[P} (31) 

where N s is the total number of corner nodes existing on the 
rotor surface. It is true, in this case, that 


= I, N,j is 
,ecto, [g] 


odes existing on the 

(?) ' f«):' 


but a subvector of the already computed global 
in Eq. (30) and is, therefore, known as this com- 



Flg. 6 Pressure integration segment over the surface of a generaity 
shaped rotor 


putational step. 

Next, consider the finite element face “s” in Fig. 6 which 
exists on the rotor surface in its displaced position. The pressure 
derivative p { J ] over this surface element can now be interpo- 
lated in terms of the known pressure derivatives at the four 
comer nodes of **s ” as follows: 


pT 


4 


f-1 



(32) 


Summation of p [ J\ as contributed by all finite elements sharing 
faces with the rotor, over the entire rotor surface yields the 
rate (with respect to “c”) at which the fluid force is exerted 
on the rotor. This can be resolved in the ct x” and “y” directions 
3F X J bF y c „ 

to produce the derivatives —— and as follows: 
be be 

\ v.i)dvdi (33) 


bFy 

be 



n y (y,Z )p,t(v>Z)G(ihZ)dvdt 


(34) 


where n x and n y are the components of the local unit vector 
that is normal to the rotor surface (Fig. 6). Also, the parameter 
GOf, £) is a Jacobian-like function for cartesian-to-local area 
transformation, and was previously derived by Baskharone 
(1979) as follows: 





by 2 /bx bz bx dz\ 2 
dZdv) + [d(dv~d v dij 


bx by 

dvb( 


bx 2 


1/2 


where the derivatives on the right-hand side are evaluated using 
the interpolation expressions (9) for a "f” value of 1.0 (Fig. 
6). Having computed the force derivatives (Eqs. (33) and (34) 
above), the requirements for computing the force coefficients 
(AT, k f C, c, Af, m) are, by reference to Eqs. (3) and (4), 
complete. 

Concluding Remarks 

The virtually deformable Finite element concept introduced 
in this paper offers a new and powerful tool with which the 
fluid/structure interaction effects can be captured. Versatility 


of the computational model was illustrated through two typical 
problem categories where a cylindrically whirling rotor was 
analyzed. In this, and all similar disturbance modes, the object 
of the perturbation analysis is the finite-element equivalent of 
the flow governing equations in the rotor-to-housing gap. 

Sample results of the current analysis are contained in a 
companion paper (Baskharone and Hensel, 1991). Under focus 
in this paper is the rotordynamic coefficient of an annular seal 
under a cylindrical type of rotor whirl. It is undoubtedly clear, 
at this point, that tilting-related, or simply moment, coeffi- 
cients of a conically-whirling rotor would equally be calculable 
using the virtual finite-element distortion concept developed 
in the current study. Utilization of the concept in this case 
would lead to the differential changes in the fluid forces along 
the rotor axis which, in turn, exert differential moments around 
two perpendicular axes at the center of the tilting motion. The 
moment coefficients are eventually computed by integrating 
these moments along the axis and relating the outcome to the 
rotor angular displacement and its time derivatives. It is further 
emphasized that the same concept is potentially applicable to 
other modes of fluid-induced vibrations. An example of these 
modes is the axial oscillation of the impeller-shroud assembly 
of the centrifugal pump to which reference was made in this 
paper. The oscillatory motion is, in this case, caused by the 
fluid forces on the front (or shroud) side of the impeller as 
well as the back face of the impeller disk. An equally detailed 
computational model that is capable of analyzing such vibra- 
tion modes is currently non-existing in the open literature. 
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A newly devised perturbation model for the fluid-induced vibration of turboma- 
chinery rotating elements is used to compute the rotordynamic coefficients of an 
annular seal First , the finite element-based solution of the flow field in the centered- 
rotor operation mode is verified and its grid dependency tested for different seal 
configurations. The rotordynamic behavior of a hydraulic seal with a clearance gap 
depth/length ratio of 0.01 , as a representative case , is then analyzed under a cylin- 
drical type of rotor whirl and several running speeds. The direct and cross-coupled 
rotordynamic coefficients dictating the rotor instability mechanism in this case are 
compared to experimental and analytical data , and the outcome is favorable. The 
numerical results are also used to discuss the validity of a common assumption in 
existing computational models in regard to the circumferential distribution of the 
perturbed flow variables in the eccentric rotor operation mode. 


Introduction 

Fluid-induced vibration of turbomachinery rotating com- 
ponents is often a result of today’s high performance demands. 
These, among others, include excessively tight clearances, high 
speeds, and highly loaded blades in the primary-flow passages. 
While these design trends intendingly maximize the overall 
efficiency, they may also lead to forced vibration of the rotating 
components. This potentially dangerous situation can also be 
the result of fluid/solid interaction mechanisms in secondary- 
flow passages such as seals, dampers and bearings. Of these, 
the problem of seal rotordynamics has been the focus of ex- 
tensive research in recent years. 

Although the primary function of seals is leakage control, 
they evidently influence the overall stability of turbomachines. 
While nonsmooth seals, of the grooved, labyrinth and honey- 
comb types (Diewald and Nordmann, 1988; Childs and Schar- 
rer, 1987; Childs, 1988) are known to produce minimum 
leakage, it is not established with certainty that they provide 
better rotordynamic characteristics when compared to smooth 
untapered seals. It is established, however, that the instability 
in all cases are dependent upon such parameters as the clearance 
depth (Childs and Scharrer, 1987), the flow preswirl (Childs, 
1988) and other geometry-related variables. The latter includes 
the taper angle, in the case of an annular seal (Nelson, 1984) 
and whether the seal is of the “look -through” or the stepped 
type (Scharrer, 1988) in the class of labyrinth seals. 

Recent advances in the seal rotordynamics area have estab- 
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lished the computational fluid dynamics tool as a strong mod- 
eling foundation. The first and foremost contribution in this 
field was that by Dietzen and Nordmann (1987), where a finite- 
difference formulation of Navier-Stokes equations was used 
in conjunction with a perturbation model to compute the ro- 
tordynamic coefficients of incompressible-flow annular seals 
(Fig. 1). The significance of that study lies in the departure it 
advocated from the bulk-flow model (Childs, 1983) in which 
details of the rotor-to-housing flow field were significantly 
simplified. A common assumption in the study by Dietzen and 
Nordmann (1987), as well as that of Childs (1983), has involved 
a single-harmonic sinusoidal variation of the flow properties 
in the tangential direction. 

It is the purpose of this paper to validate a comprehensive 
finite element-based perturbation model which the authors 
have recently devised (Baskharone and Hensel, 1991). The 
model conceptually deviates from the traditional perturbation 
analyses (Dietzen and Nordmann, 1987; Childs, 1983) in that 
the perturbation equations emerge from the flow-governing 
equations in their discrete, finite-element form as opposed to 
the differential form in all existing models. The perturbation 
in the flow field, in this case, is physically perceived as a result 
of virtual distortions in an assembly of finite elements that is 
occupying the rotor-to-housing gap. The computational ap- 
proach, as such, makes it potentially possible to analyze dif- 
ferent rotor vibration modes which may not necessarily involve 
a uniform rotor eccentricity (which leads to a cylindrical rotor 
whirl), as well as imparts the versatility and adaptability of 
the Finite-element method to this class of engineering appli- 
cations. 

In view of the substantial versatility of the current model, 
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it may appear that cylindrically whirling rotors of annular seals 
do not represent a sufficiently challenging problem, owing to 
the extensive amount of existing data that virtually covers all 
aspects of this problem. However, it is precisely the availability 
of these data that led to selecting this sample case. As dem- 
onstrated by Baskharone and Hensel (1991), the new pertur- 
bation model offers a conceptual deviation from all existing 
fluid-induced vibration models. Although the originality here 
is indeed rotordynamics-related, it is also the application of 
the finite-element technique to such excessively narrow gaps 
with a considerable magnitude of flow turbulence and inertia 
domination, that is hardly verified in the literature. Since anal- 
ysis of the zeroth-order and perturbed flow fields are both 
based on this technique, the availability of experimental data 
to verify the finite-element flow solution was a primary cri- 
terion in the selection process. As a result, the numerical in- 
vestigation in this paper asserts, among others, the finite- 
element method as a reliable means of securing the narrow- 
gap zeroth-order flow solution, with which the perturbation 
analysis is initiated. The investigation, therefore, goes beyond 
the objective of simply validating the final outcome of the 
perturbation model. 

Centered-Rotor Flow Field 

The axisymmetric flow analysis corresponding to this rotor 
position provides the zeroth-order solution in the perturbation 
model (Baskharone and Hensel, 1991). Presented here, how- 
ever, is an assessment of the stability, smoothness and accuracy 
aspects of the numerical solution as strongly tied to the level 


of nonlinearity in the flow equations, the latter being implied 
by the Reynolds number, and the manner in which the flow 
turbulence and near-wall flow structure are analyzed. One of 
two example cases under investigation is a hydraulic annular 
seal with the following dimensions and operating conditions: 

Rotor radius (a) = 31.0 mm 
Clearance ( 5 ) = 0.294 mm 
Length (L) = 200.0 mm 
Through-Flow Reynolds number (RJ= 10 4 
Tangential Reynolds number (R«.J = 2.5 x 10 4 

where the through-flow Reynolds number is based on the inlet 
through-flow velocity and the tangential Reynolds number is 
based on the circumferential velocity of the rotor surface, with 
the characteristic length being the seal clearance in both cases. 
Experimental measurements concerning this seal configuration 
were reported by Yamada (1962) in the form of a seal “re- 
sistance” parameter (X), which is defined as follows: 



where: 

A p is the static pressure drop across the seal, 
p is the fluid density, 

V z . is the seal-inlet through-flow velocity. 

Turbulence Closure and Near-Wall Flow Zone. An alge- 
braic eddy viscosity model by Baldwin and Lomax (1978) is 
used to simulate the flow turbulent behavior. According to 
this model, the effective kinematic viscosity is viewed as com- 
posed of two, molecular and eddy, components as follows: 

v e =vi + v t ( 2 ) 

In calculating the eddy component, v u the procedure assumes 
the presence of two, inner and outer, layers. In the inner layer, 
the Prandtl-van Driest formulation yields the following expres- 
sion: 

»-,,/ = / 2 lu'l (3) 

where the subscript / refers to the inner layer. The mixing 
length / in expression (3) is defined as follows: 



where: 

y is the distance normal to the nearest wall 
t w is the wall shear stress 
os' is the local vorticity 

The model switches from Van Driest formulation to that of 
the outer region at the smallest value of y for which the inner 
and outer values of the eddy kinematic viscosity are equal. The 
formulation for the outer layer is given by: 

v t,o ~ KC C p F max y mt^^KLEB (4) 

where: 


Nomenclature 


a = rotor radius 
b = housing radius 
C, c = direct and cross-coupled 
damping coefficients 
K t k = direct and cross-coupled 
stiffness coefficients 
L = seal length 
/ = mixing length 

M, m - direct and cross-coupled in- 
ertia coefficients 
p = static pressure 


Re = Reynolds number 
s = clearance gap width 
y = distance from the wall 
0 « grid refinement factor 
X = dimensionless seal resistance 
(defined in Eq. (1» 
p = viscosity coefficient 
v = kinematic viscosity coeffi- 
cient 

p = density 
r = shear stress 


fl = whirl frequency 
a) = shaft running speed 
u ' = vorticity 

Subscripts 

i - seal inlet station 
R = relative flow property 
r = radial component 
w = value at the wall 
z - axial component 
$ = tangential component 
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with y max referring to the value of y at which F miX occurs. The 
various constants in Baldwin-Lomax model are as follows: 


A* = 26, *=0.4, *=0.0168, C cp = 1.6; and CK i£fl = 0.3 

It is important to point out that in applying this turbulence 
model, which was devised for two-dimensional flow applica- 
tions, the vorticity and wall shear stress calculations were mod- 
ified in such a way that the tangential velocity component is 
taken into account. This was necessary since the flow in the 
rotor-to-housing passage is generally that of the swirling type. 

The authors experience with the foregoing turbulence closure 
during the course of this study has proven that, regardless of 
any “practical” level of mesh refinement near solid boundary 
segments, accurate values of the near-wall eddy viscosity were 
unachievable at high Reynolds numbers. The focus then was 
on the annular seal described above, for which the flow re- 
sistance, as measured by the static pressure drop across the 
seal, was persistently overestimated. Further investigation re- 
vealed that the high velocity gradients in the near-wall zone 
were improperly depicted. The situation was conclusively 
remedied by utilizing the universal law of the wall concept in 
computing the wall shear stress “r*” appearing in the expres- 
sion of the dimensionless distance y + above. 

Computation of the wall shear stress, that is associated with 
a typical node “/” (Fig. 2), is based on the near-wall zone 
treatment proposed by Benim and Zinser (1985). The as- 
sumption here is that the universal law of the wall at any wall 
location is extendible, across the wall element width, to the 
grid node (Fig. 2) at this axial location. Referring to the 
distance of node “y” from the wall by y m jn , the following 
expression for the wall shear stress is then obtained: 


y min 

«c% A P v m 


Intern — 

V. \ 1, 


for ^m in <11.6 


(5) 

for 11.6 


where: k min = -% 75, C D = 0.09, * = 0.4, £ = 9.0 

P<~D 

with K min referring to the magnitude of velocity, with the tan- 
gential component taken into account, at the interior node and 
v\ being the molecular coefficient of kinematic viscosity. Re- 
ferring to Eq. (5) above, note that the outcome of this equation 
in the case where y^ ^ 11.6 is a recursive relationship since 
t w now appears on both sides of the equation. An iterative 
procedure is executed, in this case, to compute t w . Noteworthy 
is also the fact that the node “y” (Fig. 2) in the Benim and 
Zinser model is replaced by the corner node “fc” which is 
consistent with their choice being a four-noded bilinear Finite 
element. In reality, however, neither one of these nodes would 
be sufficiently close to the wall, in view of the steep velocity 
gradient in this region. An alternate choice of this near-wall 
point was, in this case, adopted as part of the accuracy en- 
hancement process discussed next. 

Numerical implementation of the preceding turbulence clo- 
sure, including the near-wall model, is achieved with the aid 
of an array of points that is different from the primary set of 
computational nodes. Figure 2 shows an enlarged segment of 
the computational domain near a solid wall, in which the 
primary nodes in the finite-element discretization model are 
identified by hollow circles, while the points used in the eddy 
viscosity computations at the typical node are solid circles. 



Fig. 2 Refinement of the near-wall region for accurate prediction of 
the eddy viscosity 

The objective here was twofold; to estimate the cut-off location 
between the inner and outer layer with sufficient accuracy, 
and to capture the steep gradients of the flow variables near 
the solid wall. A substantial enhancement of the solution was, 
as a result, observed during the preliminary testing phase. This, 
for the major part, was due to the excellent accuracy with 
which the wall shear stress was computed, for it was generally 
possible to compute the stress at a point where the dimen- 
sionless wall coordinate (y + ) was acceptably small. 

A modification to the Baldwin and Lomax turbulence clo- 
sure was made in the current study to accommodate the angular 
motion of the rotor surface. Referring to the definition of the 
inner-layer kinematic eddy viscosity (Eq. (3)), and considering 
the case in which the layer is attached to this rotating surface, 
the vorticity and wall shear stress were both defined on the 
basis of the relative velocity (V*) as opposed to the absolute 
velocity (V), where: 

V* = V-<*7e* (7) 

where: 

u is the rotor spinning speed 

a is the rotor radius 

e e is the unit vector in the positive tangential direction 
In this case, the absolute vorticity (a>*) is expressible in terms 
of the absolute velocity (a/) as follows: 

& # = a) * — 2&J (g) 

Calculation of the eddy viscosity in the outer layer (Eq. (4)) 
as pertaining to the rotating surface was consistently performed 
using the relative flow properties. 

Introduction to the Upwinding Technique. Adoption of 
the conventional Galerkin’s method in deriving the finite-ele- 
ment equivalent of the flow governing equations was, perhaps 
expectedly, unsuccessful. This was a result of the high Reynolds 
number which caused the flow field to be inertia-dominated. 
Surprisingly, the “wiggles” appearing in the streamwise pres- 
sure distribution in this case were not entirely eliminated even 
when full upwinding of the equations of motion, in the manner 
devised by Heinrich and Zienkiewicz (1977), was implemented. 
In fact, the large magnitude of false numerical diffusion that 
was imparted to the flow field in this case also made it prac- 
tically impossible to achieve a numerical solution in the neigh- 
borhood of the experimental data. The situation was further 
worsened by the large aspect ratio of the rotor-to-housing gap 
(defined as the gap length/clearance ratio) which, in turn, 
produced high aspect-ratio finite elements. This had an adverse 
numerical effect, especially near the rotor and housing surfaces 
where the need for high resolution of the flow structure led 
to excessively narrow elements in this region. 

Corrective measures included deviation from the full up- 
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winding strategy and optimization of the elemental aspect ra- 
tio. First, in implementing the weighted-residual method 
(Zienkiewicz, 1971), terms in the error functions which result 
from the momentum equations were weighed differently in the 
process of deriving the finite element equations. In this case, 
the weight functions proposed by Heinrich and Zienkiewicz 
(1977), being applicable to the current axisymmetric flow prob- 
lem, were used in conjunction with only the convection terms 
in the momentum equations. On the other hand, the element 
shape functions were used with all other terms. This is con- 
sistent with the well-known practice in the area of finite-dif- 
ference modeling, where backward differencing is exclusively 
utilized when approximating the convection terms. Optimi- 
zation of the element aspect ratio, on the other hand, con- 
tributed to the smoothness of the flow variables as will be 
discussed in the grid dependency section of the numerical re- 
sults. 

Grid Dependency and Accuracy Assessment. Sensitivity of 
the numerical solution to the size of the finite-element discre- 
tization model, and the level of field resolution it provides, 
was first investigated. In all cases considered, the number of 
cross-flow stations between the seal inlet and exit stations was 
fixed at 48, while the number of grid lines between the inner 
and outer walls, referred to as N n was varied from 7 to 11. 
The latter family of grid lines was constructed with varying 
increments in the radial direction. To better quantify this var- 
iation, the growth of the increments from either wall to the 
mean radius was made to be of the geometric sequence type 
with a common ratio “0.” Illustrated in Fig. 3 is the geo- 
metrical effects of varying “0** from 1.0, which yields equal 
increments, to a value of 1.5 for the case where N r is fixed at 
11. Note the drastic change in the width of the wall elements 
as a result. 

Shown in Fig. 4 is a comparison between the computed seal 
resistance “X,” corresponding to different discretization pat- 
terns, and the experimental value of Yamada (1962). Exami- 
nation of this figure reveals that the accuracy of the seal 
resistance, which is largely influenced by the precision of the 
computed eddy viscosity, is generally improved by refining the 
finite-element grid near the walls, be that the result of increas- 
ing ‘TV/’ or “0.” It was in view of these results, that an 
optimum grid with N r ■ 9 and 0=1.5 was permanently fixed 
thereafter. 



housing 
3 — 2.5 
ROTOR 



1.5 



1.0 


Fig. 3 Geometrical effects of varying the refinement factor OS) on the 
finite-element discretization model 


Accuracy of the centered-rotor flow field was further in- 
vestigated, with the emphasis being on the cross-flow velocity 
profiles through the seal. Evaluation of the finite-element so- 
lution in this case was made in light of the flow measurements 
of Morrison et al. (1988) for an annular seal with a gap width/ 
length ratio of 0.034. Development of the through-flow and 
tangential velocity profiles along this seal were measured by 
Morrison et al. under a high preswirl/inlet through-flow ve- 
locity ratio of 4.2, and a Reynolds number, based on the 
clearance width and inlet velocity, of 13,280. Figure 5 shows 
the set of measured through-flow and tangential velocity pro- 
files at non-dimensional axial length ratios of 25, 50, and 75 
percent, together with the computed profiles for comparison. 
As seen in the figure, the computed velocities seem to depict 
the experimental profiles with reasonable accuracy. Note that 
the peak point of the through-flow velocity profiles, partic- 
ularly in the computed profiles, corresponds to radii which are 
higher than the mean radius due to the radial shifting of the 
flow particles as a result of the high circumferential velocity 
which, in turn, is caused by the high rotational speed of the 
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Fig. 4 Seal resistance as a function of the finlto-element model size 
and near-wall refinement 
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seal inner wall. Development of the computed tangential ve- 
locity along the seal in this figure seems to also be in close 
agreement with the flow measurements. 

Perturbation Analysis and Rotordynamie Coefficients 

The newly devised perturbation analysis (Baskharone and 
Hensel, 1991) was applied to a representative case of annular 
seal for which experimental and analytical rotordynamie data 
exist (Dietzen and Nordmann, 1987). Under focus here was 
the variation of the direct and cross-coupled stiffness, damping 
and inertia coefficients (AT, k y C, c, M , and m) of the fluid/ 
rotor interaction system. Of these, only the stiffness and damp- 
ing coefficients were reported by Dietzen and Nordmann (1987), 
and are used in the present study for verification purposes. 
The analytical data by Dietzen and Nordmann (1987) is the 
outcome of a traditional finite difference-based perturbation 
analysis in which sinusoidal variations of the field variables 
around the circumference were assumed, as opposed to cir- 
cumferential variations that are totally unrestricted in the cur- 
rent model. 

Of particular interest in assessing the new perturbation model 
were the stability-related rotordynamie coefficients of which 
experimental and analytical data were reported by Dietzen and 
Nordmann (1987) for a typical hydraulic seal, These are the 
cross-coupled stiffness and direct damping coefficients “Ar” 
and “C,” respectively, which dictate the fluid-exerted tangen- 
tial force and, subsequently, the rotor whirl. The seal under 
investigation had a clearance/length ratio of 0.0085, a rotor 
radius of 23.5 mm, and a Reynold’s number (based on the 
clearance and inlet through- flow velocity) of 4700. Operating 
speeds of 2000, 4000, and 6000 rpm were considered, together 
with inlet preswirl/through-flow velocity ratios 0.04, 0.10 and 
0.17, respectively. 

The three-dimensional grid used to analyze the distorted 
clearance gap (Baskharone and Hensel, 1991) was first optim- 
ized. In this case, the number of tangential stations “AV’ in 
this gap (Fig. 6) was varied, while maintaining the axial and 
radial station counts, as well as the near-wall refinement level, 
unchanged. Figure 6 shows the variation of the error in the 
computed force derivatives, dF x /de and 3F y /de t where ‘V’ is 
the rotor virtual eccentricity, as a result. The error in this figure 
is defined as the dimensionless difference between the force 
derivative and that corresponding to an N e value of 12, and 
the force derivatives are those associated with a zero whirl 
frequency. Based on the asymptotic shape of both curves in 
Fig. 6, the number of tangential stations was chosen to be 12. 
As a result, the RAM and CPU time consumption, regardless 
of the whirl frequency magnitude, were roughly 1220 K. bytes 
and 12 minutes, respectively, on the SX2 supercomputer, per 
each whirl frequency. 

Shown in Fig. 7 is a comparison of the computed rotor- 
dynamic coefficients with those reported by Dietzen and Nord- 
mann (1987). The latter data set contained measurements of 
the direct and cross-coupled stiffness coefficients, but only the 
direct damping coefficient. Also contained in this set were the 
numerical results obtained by Dietzen and Nordmann, for all 
four of the stiffness and damping coefficients, using a tradi- 
tional perturbation analysis in which the finite-difference 
method was utilized. As seen in this figure, the computed 
results are generally in close agreement with the experimental 
and analytical data as far as the cross-coupled stiffness and 
direct damping coefficients are concerned. These two coeffi- 
cients, as discussed earlier, are representative of the fluid- 
induced tangential force perturbation and, as such, dominate 
the mechanism of rotor instability. The cross-coupled stiffness 
coefficient “At” is, by reference to Fig. 7, positive, therefore 
destabilizing, and increases with the rotor speed. The direct 
damping coefficient "C,” on the other hand, provides a sta- 
bilizing effect which, according to the computed results, in- 
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Fig. 6 Sensitivity of ths force derivatives to the number of grid planes 
around the circumference 
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Fig. 7 Comparison of the direct end cross*coupled rotordynamie coaf 
ficients with experimental and analytical data 


creases slightly with the operating speed. Examination of Fig. 
7 reveals, however, that the predicted values of the direct 
stiffness coefficient %t K ” and cross-coupled damping coeffi- 
cient “c” are significantly different from those computed by 
Dietzen and Nordmann (1987). This, for the major part, is 
caused by differences in the unperturbed flow fields including, 
in particular, the radial profiles of the flow properties along 
the seal axis as explained next. 

Evolution of the direct stiffness and cross-coupled damping 
coefficients in Fig. 7 can be traced back (Baskharone and 
Hensel, 1991) to a single source, namely the radial component 
of the fluid reaction force on the rotor. As reported by Hensel 
(1990) for the same seal, this force component is highly sensitive 
to the inlet through-flow velocity profile including, in partic- 
ular, the inlet boundary layer thickness. For instance, Hensel 
reported an average reduction of nearly 21 percent in the direct 
stiffness coefficient “Af ” as a result of arbitrarily reducing the 
inlet boundary layer thickness, on both the rotor and housing 
surfaces, from 15 percent of the clearance width to 7.5 percent, 
which is the current value, with the cross-coupled stiffness and 
direct damping coefficients remaining practically unchanged. 
Such sensitivity would naturally be anticipated in seals, such 
as the present, where the axial length is too small to permit a 
fully developed flow even at the exit station. It is therefore 
true that the computed values of li K ” in Fig. 7 are primarily 
dependent on the inlet boundary layer thickness assumption 
which, in turn, was necessary since no measured thickness was 
available. On the other hand, there is no clear evidence that 
this assumption had any significant effect on the trend of “A"” 
with the rotor speed in Fig. 7. Nevertheless, it should be noted, 
in reference to this figure, that the experimentally-determined 
peak of 14 AT” at a rotor speed of 4000 rpm is as absent in the 
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results of Dietzen and Nordmann (1987) as it is in the current 
results despite the major difference in simulating the flow 
turbulence in both studies. 

Versatility Assessment of the Current Perturbation 
Model. Thorough examination of the current numerical re- 
sults for the simple annular seal under consideration appeared, 
as anticipated, to confirm the validity of the single-wave sin- 
usoidal forms of flow perturbations around the circumference. 
The flow variable under focus was the rotor-surface pressure 
perturbation (dp/de) at the middle seal section which, through 
use of the computed nodal pressures, was cast in the following 
Fourier series form: 

= <*o + 2 cos ("0 ) + 2 b » sin ("0) ( 9 ) 

n n 

where the subscript **/” refers to the inner (or rotor) radius 
and “0” is the angular coordinate. The number of sine and 
cosine waves * V* in expression (9) was arbitrarily set equal to 
five. The ratios (a n /a x ) and (b„/b\) were then computed. All 
of these ratios came to be less than 0.5 percent for 2 < n < 
5 regardless of the whirl frequency/spinning speed ratio. This 
simply means that the common form of field variable pertur- 
bation, namely: 

(P = (P c cos 0 + ( P s sin 0 (10) 

where (P stands for the perturbation of any flow variable, is 
seemingly as accurate as can be. This functional form is re- 
cognizably a built-in assumption in any existing perturbation 
model in this rotordynamics area (e.g., Dietzen and Nord- 
mann, 1987; Childs, 1983). 

Complex clearance gaps and/or vibration modes, however, 
are where approximations of the type indicated above may 
become inapplicable. Consider, for instance, the problem of 
a shrouded pump impeller (e.g., Baskharone and Hensel, 1989) 
where the rotor, being the shrouded impeller in this case, under- 
goes a conical type of whirl as a result of an angular eccentricity 
of the rotor axis. As indicated by the authors (Baskharone and 
Hensel, 1991), it is this problem configuration to which the 
virtual finite-element distortion concept, being the foundation 
of the current perturbation model, is uniquely applicable. In 
this case, deformations of the finite elements, due to the an- 
gular rotor eccentricity, will give rise to differential changes 
in the rotor surface pressure which, in the end, produce fluid- 
exerted moments M x and M y around the center of conical whirl. 
Future application of the current analysis to this type of rotor 
eccentricity may indeed make it possible to assess the validity 
of the simple sinusoidal distribution (indicated above) under 
such circumstances. The authors suspect that such simple dis- 
tribution would probably be unrealistic in the tooth-to-tooth 
gap of the labyrinth seal which may very well be part of the 
leakage passage (e.g., Baskharone and Hensel, 1989), as well 
as other subdomains where flow separation and recirculation 
occur. Stated differently, it is in critical subdomains, such as 
these, that the coefficients “a*” and “b n ” in expression (9) 
may not be necessarily zero for n > 1 as is generally assumed 
to be the case in virtually all existing perturbation models of 
fluid-rotor interaction. 

Concluding Remarks 

The intention in this paper was to verify a new and cate- 
gorically different methodology in the area of fluid-induced 
vibration, as well as validate its fluid mechanics foundation. 
Both tasks were successfully accomplished using representative 


annular seal configurations. The direct and cross-coupled ro- 
tordynamic coefficients, which comprise the ultimate goal of 
the study, were in close agreement with available experimental 
and analytical data, particularly those which are directly linked 
to the rotor instability. While the strategy in the current per- 
turbation model was clearly different from all existing per- 
turbation models in this field, an effort to assess the versatility 
of the new model, by comparison, was made and the funda- 
mental differences identified. Since limitations on the rotor- 
to-housing clearance gap in the model are practically non- 
existing, it is presently the plan to apply the model to shrouded 
impellers of the Space Shuttle Main Engine turbopumps for 
which the destabilizing forces are developed in the shroud-to- 
housing secondary-flow gap. Release of the computer code, 
as part of an upcoming NASA contractor report, will take 
place upon completion of this final research phase. 
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MOMENT COEFFICIENTS OF INCOMPRESSIBLE-FLOW SEALS 
WITH CONICALLY WHIRLING ROTORS 

E. A. Baskharone and S. J. Hensel* 

Mechanical Engineering Department, Texas A&M University, College Station, TX 77843, U.S.A. 

( Received 1 May 1990, and in revised form 13 September 1990) 

Abstract — A computational model is -developed to investigate the dynamic behavior of a rotor that 
is in contact with a fluid, to infinitesimally small disturbance modes. The rotor here is that of an 
annular seal, commonly used in turbomachinery applications, with an incompressible flow in the 
annular clearance gap. Under consideration is an angular excitation of the rotor axis, coupled with a 
whirling motion around the housing centerline at a finite whirl frequency. The fluid response in this 
case is quantified in the form of direct and cross-coupled moment coefficients, which constitute a 
measure of the stiffness, damping and inertia effects of the rotor/fluid interaction. Uniqueness of the 
computational model stems from the manner in which the rotor eccentricity is physically perceived 
and subsequently incorporated. It is first established that the fluid reaction components are the 
result of infinitesimally small deformations of varied magnitudes that are experienced by an 
assembly of finite elements in the rotor-to-housing gap as the gap becomes distorted due to the rotor 
virtual eccentricity. The idea is then cast into a perturbation model in which the perturbation 
equations emerge from the flow-governing equations in their discrete finite-element form, as 
opposed to the differential form, which is traditionally the case. The computational model is 
potentially applicable to a wide range of rotor-to-housing pasage configurations where ellipticity of 
the passage flow field is too significant to ignore. The numerical results are compared to those 
obtained through an existing bulk-flow model which is particularly suited for the leakage passage in 
annular seals and other similarly simple passage configurations. 


1. INTRODUCTION 

Seals are conceptually leakage control devices which are extensively used in turbomachines 
to minimize the secondary flow rate in clearance gaps. Of these the straight type (Fig. 1) is 




•) CYLINDRICAL WHIRL 



b) CONICAL WHIRL 


Fig. 1. Cylindrical and conical whirl of an annular seal rotor. 
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often utilized in shrouded-impeller turbopumps to reduce leakage over the impeller front 
face. However, it has long been recognized (e.g. Black [1]) that the seal rotordynamic 
characteristics also have a primary influence on the stability of the entire turbomachine and 
are, therefore, crucial in the life estimation process. This motivated several researchers to 
compute [2-5] and measure [5, 6] the stiffness, damping and inertia coefficients resulting 
from the fluid/rotor interaction. Despite the difference in methodology among com- 
putational analysts in this area, the definite majority of numerical models existing today are 
those concerned with the force coefficients associated with cylindrical whirl of the rotor axis 
(Fig. la) ' 

The study by Childs [3] was the first and remains the only analytical attempt to quantify 
the moment-related coefficients in annular seals. This was a perturbation analysis based on 
the Hirs’ bulk-flow model [7] where the rotor tilting motion was coupled with the rotor 
whirl around the centered position (Fig. lb). Childs’ conclusions were later validated by 
Kanemori and Iwatsubo [6] for a long seal where the fluid-induced moments would 
naturally be of predominant influence on the seal dynamic behavior. The experimental 
study by Kanemori and Iwatsubo was, nevertheless, incomplete for it was only concerned 
with the fluid-exerted moment resulting from a uniform lateral eccentricity (cylindrical 
whirl). Figure 2 illustrates this type of rotor excitation which, considering the axial gradient 
of the rotor surface pressure, indeed leads to tilting moments, especially in long seals. The 
impact of these moments becomes negligible as the seal length is made smaller, a case that is 
best represented by the so-called neck ring seal. Figure 3 shows this seal as part of the 
secondary flow passage in a pump stage that is similar to the booster impeller stage of the 
Space Shuttle Main Engine oxidizer turbopump. 

The current computational study deals with the direct effects of an angular eccentricity 
that is simultaneously coupled with a whirling motion of the rotor at a finite whirl frequency 
Q as depicted in Fig. 4. The analysis is based on, and intendingly advocates, what may be 
referred to as the virtually distorted finite element concept. Aimed at the direct and cross- 
coupled moment coefficients of the seal, the computational model is particularly applicable 
to short seals with length/diameter ratios that are generally less than unity. 



Fig. 2. Fluid-exerted moments in long seals as a result of the pressure differential across the rotor. 
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Fig. 3. Neck ring seal as part of the secondary passage in shrouded-impeller pumps. 


V 




¥ 



Fig. 4. Fluid-induced moments on a conically whirling rotor and definition of the coordinate axes. 


2. ANALYSIS 

Components of fluid/rotor interaction 

Figure 4 shows a whirling rotor which is experiencing an infinitesimally small angular 
eccentricity 9 around the center of tilting motion. Referring to the whirl frequency by ft, the 
eccentricity can generally be resolved at any given time t into two angular displacement 


v 
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components around the x and y axes as follows: 

9 X = -cos(fir), 0., = - sin (fit) (1) 

z z 

with z being the distance from the center of tilting motion (Fig. 4). The fluid reaction, in this 
case, is in the form of differential moments <5Af , and 5M y around the x and y axes, for which 
the following linearized relationship [2-5] applies: 



where 6 X and 9 y are the components of angular displacements of the rotor axis around the x 
and y axes, respectively. The direct (. K , C f A?) and cross-coupled (£ c, m) rotordynamic 
coefficients, in the square matrices above, are the seal stiffness, damping and inertia 
coefficients under an angular rotor excitation. The differential moments 6M X and SM y in 
equation (2) represent the fluid reaction to this excitation, where: 


5M X — J zSF y \ and 

(3) 

SM y = f zSF x 

(4) 

with SF X and 5F y being the differential changes in the forces exerted on the rotor due to the 
infinitesimally small shift in pressure on the rotor surface as a result of the angular 
eccentricity (Fig. 4). With no lack of generality, consider the special location of the displaced 
rotor axis in a plane that is at a distance z from the center of the tilting motion as shown in 
Fig. 1(b). Using equation (1) in this case yields the following expressions for the angular 
location, velocity and acceleration of the rotor axis: 

0 i = i = 6> ( 0 y = O 

(5) 

0 i = o, 0 = - = £10 

Z 

(6) 

Q x = e -= — fl 2 0, 0, = O. 

(7) 

Use of expressions (5)— (7) in the matrix equation (2) gives rise 
equations: 

to the following two 

-*-*+«■* 

(8) 


(9) 


Determination of the stiffness, damping and inertia coefficients in equations (8) and (9) 
requires computation of dM x /86 and 8M y /89 at three or more different values of the whirl 
frequency fl Interpolation of these two derivatives as quadratic functions of D, using curve 
fitting techniques, leads to the evaluation of etc. by simply equating corres- 

ponding terms on both sides of equations (8) and (9). 

The procedure to determine 8MJ88 and 8M y /86 above is initiated by securing a finite- 
element solution of the undisturbed flow field in the rotor-to-housing annular gap for the 
centered-rotor operation mode. Next, a virtual angular eccentricity 9 (Fig. lb) is introduced. 
As a result, the finite-element nodal points undergo virtual lateral displacements of different 
magnitudes depending on the spatial location of each node (Fig. lb). Analysis of the 
perturbed set of finite-element equations under such virtual distortion yields, among other 
variables, the rate at which the rotor surface pressure varies with the angular eccentricity 
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(i.e. dp/dd). The coefficients dM x /d8 and dM y /dd are, in the end, determined by integration 
along the rotor axis. 

Zeroth-order flow field 

The undisplaced rotor flow problem is defined in a frame of reference that is attached to 
the rotor (Fig. 4) and is whirling around the housing centerline at the rate of fi. It may 
appear at this point that selection of these rotating axes serves no purpose other than 
complicating the flow-governing equations. However the reader is reminded that this 
selection leads to a steady relative flow field (as viewed in the rotating frame of reference) 
once the rotor displacement occurs. Prior to this displacement, the coordinate axes are of 
the simply rotating type (Fig. 5) and the relative flow field (referred to as the zeroth-order 
flow field) is also steady. 

Definition of the coordinate axes in the manner described above, is not intended to imply 
the need to solve the flow-governing equations in the physically distorted clearance gap (e.g. 
Dietzen and Nordmann [4]). The fact, however, is that this choice, despite the complex and 
three-dimensional flow equations it produces (as will be seen next), makes it systematic to 
develop the perturbation model as it ensures consistency of the unperturbed and perturbed 
flow fields, of which only the earlier is physically axisymmetric. Note that this consistency 
requirement can alternately be ensured by solving the simple axisymmetric flow field in the 
undistorted passage, and then casting this “zeroth-order” solution in the frame of reference 
defined above [8]. This procedure was actually implemented during execution of the 
computer code for the substantial reduction in the core and CPU time consumption it 
provides. 

Viewed in the preceding frame of reference, the swirling flow in the rotor-to-housing gap 
is assumed adiabatic, incompressible and generally turbulent. The momentum and mass 
conservation equation can therefore be expressed in the rotating frame of reference (noting 
that Cl in Fig. 5 is negative) as follows: 


.du du *du ^ _ 1 dp 

&T- + + w— + 2Clv - Cl 2 x = + v,V 2 u 

ox dy dz pdx 


+ 2 


dv { du dv { fdu dv \ dv t ( du dw 

dy + dx ) + dz \dz + dx 


dx dx 


+ 


( 10 ) 


.dv dv A dv _ . 1 dp 

&T- + + W— - IClu - Cl 2 y = 4- v e V 2 p 

dz pdy 


dx dy 
+ 


^dv { dv dv t fdu dv t /dv <5w\ 

dy dy dx + dx ) + dz + dy ) 


dw dw 


. dw 


1 dp dv t dw dv t /du dw 

+ C— -b w— = -f + v e V 2 w + 2-^— + w 

dx dy dz pdz dz dz 


dv t fdu dw \ 
dx \5z + dx ) 


(ID 

( 12 ) 



Fig. 5. Definition of the rotating axes in the undisturbed-rotor operation mode. 
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du dv dw 
dx + dy 3z 


= 0 , 


( 13 ) 


where u, i\ w are the relative velocity components, p is the static pressure, Q is the rotational 
frequency of the coordinate axes (to be set equal to the whirl frequency as the rotor enters its 
disturbed operation mode) and v„ v„ v e are the molecular, eddy and effective components of 
kinematic viscosity coefficients, respectively, where v e = v, + v t , with the symbol O de- 
signating values that are carried over from a previous iteration or an initial guess for the 
purpose of successively linearizing the momentum equations. As seen, the whirl frequency O 
is part of these equations, which is a result of including the centrifugal and Coriolis 
acceleration effects produced by the axes rotational motion. 

Simulation of the flow turbulence is based on the algebraic eddy-viscosity turbulence 
model by Baldwin and Lomax [9]. However, an adjustment in this vorticity-based model 
was implemented whereby the relative, as opposed to the absolute, vorticity was used to 
calculate the mixing length in the fluid layer that is adjacent to the rotor surface. Equally 
important in implementing the model was the analysis of the near-wall zone which was 
conceptually based on the Benim and Zinser approach [10]. 

Boundary conditions 

Referring to the computational domain in Fig. 6, appropriate boundary conditions are 
needed over all boundary segments as a result of uncompromising the flow ellipticity in the 
current flow model. These include known profiles of inlet velocity components and zero 
streamwise diffusion of the exit velocity vector. The latter replaces a zero velocity gradient 
exit condition which would imply a fully developed flow that may or may not prevail at the 
exit station depending on the seal length. The last category of boundary conditions 
concerns the rotor and housing surfaces, and is created by the rotation of the coordinate 
axes. Referring to Fig. 5, the housing surface, to an observer in the rotating frame of 
reference, will no longer appear stationary, but will rather possess a relative velocity 
component Qb as indicated in the figure, where b is the housing radius. Considering the case 
of what would amount to forward whirl (as the rotor axis becomes eccentric), the rotor 
surface will appear to the same observer as rotating at an angular speed that is less than the 
rotor operating speed co by the amount £1 The rotor surface velocity in this case is (w - fi)a, 
where a is the rotor radius. 


Finite-element formulation 

Numerical solution of the centered-rotor flow field is based on the Petrov-Galerkin 
weighted-residual principle in which the weight functions are so chosen to incorporate the 
upwinding concept. An effective tool in alleviating spurious pressure modes, this concept 
has been advocated by many analysts in handling highly convective flows, such as the 
present. The concept is equally opposed by others due to the false “numerical” diffusion it 
imparts to the finite element equations. The desire to minimize this false diffusion in the 
current flow model led to selecting a mixture of the weight functions devised by Heinrich 
and Zienkiewicz [11] and those used in the conventional Galerkin method in such a way 
that upwinding is applied only to the nonlinear convection terms on the left-hand side of the 
momentum equations (10)— (12). 



Fig. 6. Finite-element discretization model (drawing not to scale). 
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The finite-element model is constructed with the 20-node quadratic element in Fig. 7 as 
the discretization unit. The procedure is initiated by defining the element shape and the 
Cartesian-to-local system mapping relations. Within a typical element (Fig. 7), let the 
spatial coordinate conversion be defined, in terms of the nodal values, as follows: 


20 20 20 

* = Z N t x i' y = Z N i z = Z 

i= i ;= i i = l 


(14) 


where N t are quadratic “shape” functions [12] of the local coordinates (, t] and 
Conversion of the spatial derivatives is, in this case, defined as follows: 


(15) 


( o 


d ) 

dx 


dt 

d_ 

dy 

> = in < 

-> 

dr, 

d 


d 

< dz J 


Ik) 


where the components of the matrix [T] are as follows: 

^ 1 ( dN, 4° dN ( |2 dN, ™ dN ( \ 

T “ - ITT ( Tt z ‘ - w y ‘ ) 

_ 1 / 8 N r 20 BN, 20 BN, 20 8N . \ 

12 \J\{& bz y, & Zi dt; dt, Z, J 

13 |J| ( 8C y ‘& dr, Z ‘ dr, y, ,ii dt, Z, J 

21 |J| ( X 'i?i dr, Zi dr, d£ Z 7 

_ 1 ( 20 5W, 20 aW| 20 20 \ 

722 - mV .5 T x ' & IT ■ *' " £ af x, ,? 1 "ar z, 7 

23 |J| ( ,5, dr, Xi ( ?! ac Z ' ac X ',?, dr, Z 7 

„ l { ™ 8 N, ™ dN, 20 20 3 jv. \ 

31 I y I V ,5-1 dr, Xi ,5 1 5$ ,5, dr, y ‘) 

„ 1 / ™ BN, ™ BN, 20 20 ajV( \ 

732 ■ \T\ ( £ If x ‘ ir J y ‘ ■ & af X| & * j 

_ 1 / 20 dN, 20 5Af| 20 20 8Nj \ 

733 = \T\ V & IT x ' & aT ■ y< " £ aT x ' & y, J 

with | J | being the Jacobian of transformation. 

Next, the flow variables are interpolated within the typical element in a similar fashion. 
Guided by the Ladyshenskaya-Babuska-Brezzi compatibility requirements [13], as ap- 
plied to the current problem, the velocity components and pressure are expressed as follows: 


20 20 20 8 
u = Z N i u « v = Z N t v t> w = Z N i w i* p = Z M rPk . 

i -1 1 i -1 *«1 


(16) 


where the interpolants M k appearing in the pressure expression are the linear shape 
functions associated with the finite element corner nodes [12]. 

The error functions produced by the unperturbed flow governing equations as a result of 
the interpolation expressions above, are then made orthogonal to a special set of weight 
functions throughout the finite element. In constructing these functions, the so-called 
error consistency criterion of Hood and Taylor [14] is implemented, whereby the element 
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• VELOCITY IS A DEGREE OF FREEDOM 
® VELOCITY AND PRESSURE ARE DEGREES OF FREEDOM 



Fig. 7. Quadrilateral isoparametric finite element used for field discretization. 


linear shape functions, M h are used in conjunction with the continuity equation. The weight 
functions used with the momentum equations, on the other hand, are of the upwind form W { 
described by Heinrich and Zienkiewicz [11] when applied to the convection terms, and are 
identical to the element quadratic shape functions N t for all other terms. 

The orthogonality conditions, stated above, constitute the elemental set of eqations, 
which can be written as follows: 

Oi]{«} + MM + MM} + M{p} = {*i} (17) 

M{“} + MM + MM + MM = (M (18) 

MM + [a 10 ]{t>} + [a u ]{w} + 0 12 ]{p} = {b 3 } (19) 

[ a l 3 ]M + [ a l4]M + [ a 15]{ w } = 0> (20) 

where 

«i)w = I { WFmF(Nj) + G(Nj)G(Nj) + H(N,)H(N j )] 

J V<" 

- ^[G(* t )G(Ar,) + H(A)H (Nj) + 2F(A)F(N,)] 

+ + CG(Nj) + wH(Nj)] }dV 

a 2 )ij = I - N ( [G(* t )F(ty) - 2nNj-\dV 

J »*> 

a 3 )ij = f - ^H(A)F(N,)dF 

J F*> 

a A \,k = f -N t F(M k )dV 

J 

a s)i,j = - N t {F(i t )G(Nj) + 2ClNj-\dV 

a 6 )i.j = |^{v t [FW)F(^) + G(N,)G(Nj) + H^jH^)] 

- ^[F(*,)F (Nj) + H(0,)H (Nj) + 2G(A)G(N,)] 

+ W t [tiF(Nj) + 6G(Nj) + WH(N ; )]}dF 
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where i = 1, 2, . . . , 20; j = 1, 2, . . . , 20; k = 1, 2, . . . , 8; the superscript ( A ) refers to a 
value that is known from a previous iteration or an initial guess, d S is the differential 
element of surface area (Fig. 7) and n is the local unit vector normal to the boundary. Also, 
the volume differential d F and the operators F, G and H are defined as follows: 


dV= |J|dCd»fd(J 


8 

8 8 

(21) 

F -r n ^+7’ 12 - + r 13 - 

8 

8 8 


G -F Jl ^+r 22 -+r 23 - 

(22) 

»- t 4 +t 4 +t 4- 

(23) 


It should be pointed out that expanding the finit element equations in terms of the operators 
F, G and H, defined above, is consistent with, and indeed simplifies, the procedure leading 
to the perturbed version of these equations, as will be seen later in this section. 

Equations (17H20) can be rewritten in the following compact form: 

MW = {*>}• (24) 

where the vector {<£} contains the nodal values of velocity components and static pressure 
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that are associated with the typical finite element. The global form of the elemental set of 
equations (24) is obtained by assembling all finite element contributions, and introducing 
the various boundary conditions, at which point the equations are expressible in the 
following form: 

LOW « {£}. (25) 

Perturbation model 

Referring to equations (8) and (9), determination of the rotordynamic coefficients reduces 
to the fundamental problem of computing the rate at which fluid-exerted moments in the 
clearance gap are developed with respect to an infinitesimally small angular eccentricity of 
the rotor axis. This is the problem under focus in this section. 

The current model is centered around the manner in which the finite-element equations 
(25) are altered as a result of the rotor virtual angular eccentricity 6 which leads to the rotor 
conical whirl depicted in Fig. 1 b. Referring to this figure, this infinitesimally small eccentri- 
city will create a linearly varying displacement, e, of the rotor axis where: 

e = 6z (26) 

with z being the distance from the center of the tilting motion. Under such lateral 
displacement, each finite element will, as a result of virtual distortion, yield a perturbed set 
of equations. Let the assembled form of these equations, by reference to (25), be as follows: 

(01] + 004] MW + {<W) = ({^} + 0{B}). (27) 

Subtracting (25) from (27), the following expression for <W> is obtained: 

(W = 0[/l]- 1 ({B}-[>4]{<i>» (28) 

or, 

d{<D} (M>\ , _ 

■V'lsbv - ™' l29) 

Equation (29) implies that two new arrays, namely the matrix [/l] and the vector {B}, are 
needed to compute the rate at which the flow variables (namely the velocity components 
and pressure) vary with respect to the rotor angular eccentricity 6. These arrays are 
primarily a result of finite elements distortion effects that are associated with the rotor 
displaced position (Fig. 8). The arrays [1] and {B} are derived next as part of the procedure 
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Fig. 8. Destruction of the flow axisymmetry as a result of the rotor eccentricity. 
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to compute the rotordynamic coefficients. To do this, a typical distorted element will first be 
analysed, and the elemental contributions [a] and {F) derived. The global arrays [A] and 
{£} can subsequently be obtained by assembling [a] and {F} among all finite elements. 

Analysis of a distorted finite element in the whirling frame of reference. Consider a 
distorted finite element in the clearance gap that corresponds to the displaced rotor 
operation mode (Fig. 8). A typical node i of this element will now be displaced by an amount 
that is a function of the local value of the rotor lateral displacement (which is a function of 
the axial coordinate) and the node location in the cross-flow plane. Referred to the axes in 
their displaced location, the new nodal coordinates are: 

X, = x„ Y i ^y i + k i 6 and Z f = z,. (30) 

The quantity A, attains its maximum value on the housing (where it is equal to - z t ) and 
vanishes for nodes on the rotor. The reason is that the housing, to an observer in the 
whirling frame of reference, is the surface that undergoes the virtual displacement (note the 
positive direction of the y-axis in Fig. 8). The local-to-Cartesian transformation Jacobian 
|J| can correspondingly be written for a distorted finite element as 1.7,1 where: 

|d,| = |J| + 0|/| (31) 

in which the Jacobian |J| is that of the undistorted element, and |/| is as follows: 


|/| = 


A new spatial derivative operator, that is equivalent to that in equation (10) can now be 
derived as follows: 
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(32) 


where the matrix [F] is as described in equation (15) while the matrix [P] is defined as 
follows: 
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Substituting expressions (31) and (32) into the flow equations (10)— (13) and re-applying 
Galerkin’s weighted-residual procedure, the distorted element equations can be written 
(upon separation of terms containing 6) as follows: 


(M + 0[fl])( w + W}) = m + e{b}). 


(33) 


where the matrix [a] and vector {h} result from the distortion in the element shape. The 
vector {<50} in this case contains the differential changes in the nodal values of velocity 
components and pressure that are created by the virtual eccentricity. 

Noting that the matrix [d] is similar in construction to [a] in equation (24), the former 
can consistently be viewed as composed of submatrices [fliHflu], which correspond to 
those in equations ( 1 7)— (20), and are defined as follows: 
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In these expressions, a variable with an asterisk has the same form as that in equations 
(17)— (20), with the exception that the Jacobian | J\ is now replaced by |7| as defined in 
equation (31), also dV = | JJdCd^d^. The vector {b} in equation (33) can similarly be 
defined by the subvectors {Z^}, {b 2 } and {fc 3 }, in consistency with equations (17H20), as 
follows: 


n/x N,x t )\J\di;dr,dS 
b 2 )i = n 2 f N,( X N, yi ) |J|dCdifd« 

J^.l \l-l / 


-°i, 


NA £ N t X t )dV 

y«) \f« 1 / 


fall = 0, 


where i and j vary from 1-20, while k varies from 1-8. With the operators F, G and H being 
those defined by equations (21)— (23), the new operators I, J and K are defined as follows: 

l=Pii M +Pii ^ +Pi3 ^ (34) 

J = P2i ^ +/>22 ^ + p “^ (35) 

K = p "ft + p “fy + p ”jz (36) 

with the matrix [P] being that defined earlier in this section. 

Calculation of the rotordynamic coefficients. Equation (29) can be rewritten in the 
following detailed form: 

£-ss W) -w-'mi-mm < 37 > 

IwJ 

where {(/}, {V}, {W} and {P} are vectors containing the nodal values of the velocity 
components and pressure, respectively, with the overbars signifying quantities that are 
associated with the perturbation in the rotor-to-housing flow field. Of all the pressure nodal 
values in the vector {P} above, define a subvector {& } as composed of the pressure values 
at the rotor surface nodes, i.e. 

{&} = {p t ,i=l > N,}<= {P} (38) 
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where N, is the total number of corner nodes existing on the rotor surface. It is true, in this 
case, that {dlP/dd} = {(8p/86)j, i = 1 , } is but a subvector of the already computed global 

vector {5<J >/30} in equation (37) and is, therefore, known at this computational step. 

Next, consider a finite element face “s” that exists on the rotor surface (Fig. 9). The 
pressure derivative p[l’ over this surface element can now be interpolated in terms of the 
known pressure derivatives at the four corner nodes as follows: 


P?= I if, 0 

i = j 


ffy) 

{86 y, ' 


(39) 


Once integrated over the proper differential segment of the rotor surface, the function p_ t 
(above) yields a differential force which in turn, can be resolved in the x and y directions. 
Taking the moments of these forces around the center of tilting motion in Fig. 2 and 
summing the result over all finite elements sharing faces with the rotor surface, the rate (with 
respect to 6 ) at which the fluid reaction moments are exerted can be expressed as follows: 


SM X 

86 

8M y 

86 



Op, tin, OG(n, £)d>/d£ 
*»*(>/, Op, e (n, OG(n, OdndO 


(40) 

(41) 


where n x and n y are the components of the local unit vector that is normal to the rotor 
surface. Also, the parameter G{t], <!;) is a Jacobian-like function for Cartesian-to-local area 
transformation, and was previously derived by Baskharone [15] as follows: 

q{„ = \ ( 8y8z 8y8z \ 2 / 8x 8z 5xdz\ 2 / dxdy dxdy\ 2 ~| 1/2 

L \8n8Z 8Z8ti) \dzty~fridz) + \fy8Z~8Z8n) J ’ 

where the derivatives on the right-hand side are evaluated using the interpolation ex- 
pressions (16) for a C value of 1.0 (Fig. 9). 

At this point, the requirements for computing the moment coefficients (K, (c, C, c, M, m) 
are, by reference to equations (8) and (9), complete. The process, as described earlier, consists 
of computing the two moment derivatives in equations (40) and (41) at three or more whirl 
frequencies fi to create a quadratic function of 8MJ86 and 8M y /86 in ft. Substitution in 
equations (8) and (9) finally leads to the moment coefficients K, k, C, . . . by inspection. 


3. RESULTS AND DISCUSSION 

A hydraulic annular seal with a clearance/length ratio of 0.0085 was selected for the 
purpose of verifying the current computational model. The force coefficients of this seal 
were determined by Dietzen and Nordmann [2], and later verified by Hensel [8]. This seal 
has a rotor radius of 23.5 mm and a Reynolds number (based on the clearance and inlet 
through-flow velocity) of 4700. Operating speeds of 2000, 4000 and 6000 rpm were 



Fig. 9. Integration of the pressure forces on the rotor surface. 



Moment coefficients of incompressible-flow seals 


165 


considered in this case, together with preswirl/inlet velocity ratios of 0.04, 0.10 and 0.17, 
respectively. 

Shown in Fig. 10 is an assessment of the zeroth-order undisplaced-rotor flow solution 
under the operating conditions cited above. Compared in this figure is the computed 
nondimensional pressure drop across the seal, to that obtained through Yamada’s experi- 
mental correlation [16]. This correlation is probably the most accurate empirical means of 
determining the seal resistance in terms of the seal geometry, rotor speed and preswirl ratio. 
As seen in Fig. 10, the computed pressure drop is in good agreement with Yamada’s 
predictions as the average deviation between the two data sets in the seal operation range is 
approximately 5%. 

The computed rotordynamic coefficients are compared to those obtained through Childs’ 
bulk flow analysis [3] for the same seal configuration and operating conditions. Worth 
noting here is the fact that Childs’ simplified model has reportedly been successful [6] when 
applied to simple clearance gaps where zones of flow separation and recirculation, such as 
those in labyrinth seals, are non-existing. Absence of these complications provides an 
environment where Childs’ idealization of the clearance gap flow, including parabolized 
flow equations and unsimulated shear stresses away from solid walls, would have a 
minimum impact on the numerical solution. Note that the current computational model is 
conceptually capable of handling such flow complications since the model is based on a 
totally elliptic flow field in the clearance gap. However, it is felt that the algebraic turbulence 
closure in the current analysis may not be appropriate for seal configurations which are 
inherently associated with massive flow separation and recirculation (e.g. the labyrinth seal 
category). Under such circumstances, a two-equation turbulence closure (e.g. Benim and 
Zinser [10]) would provide a more accurate means of modeling the flow field. 

Figure 11 shows a comparison between the computed moment coefficients and those 
obtained through Childs’ bulk-flow model [3]. These are the direct ( K , C) and cross- 
coupled (C, c) stiffness and damping coefficients defined in equation (2). The direct inertia 
coefficient (A?) in this equation was extremely small in the two sets of results, making the 
magnitudes of this coefficient comparable to the numerical errors in both computer codes. 
On the other hand, it is known [6] that straight annular seals, such as the present, are 
categorically incapable of displaying any cross-coupled inertia coefficient (m) whatsoever. 


— 1 # # # CURRENT ANALYSIS 

YAMADA’S EXPERIMENTAL CORRELATION^] 
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Fig. 10. Comparison of the computed seal pressure drop with that obtained through Yamada’s 

experimental correlation formula. 
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WWW CURRENT ANALYSIS 
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Fig. 11. Comparison of the computed moment coefficients with those computed through Childs’ 

bulk flow model. 


Interpretation of the numerical results in Fig. 1 1 is made by reference to the rotor 
components of motion illustrated in Fig. 1 and the moment derivatives in equation (9). It is 
seen in the figure that conical whirl, which is the instability mechanism under consideration, 
is created and sustained by the moment component M y , and that a positive value of dM y /S6 
indeed creates an aggravating, as opposed to restoring, effect in the case of forward whirl 
(where the rotor whirl is in the same direction as the spinning speed c u), which is the case 
depicted in Fig. lb. With this in mind, equation (9) suggests that such destabilizing effect 
would be created by positive and large magnitudes of £ and would diminish at sufficiently 
large and positive magnitudes of C. Examination of the results in Fig. 11, in light of this 
discussion, reveals that the seal under consideration becomes more stable as the rotor speed 
co is increased. This seems to consistently be the case in the speed range under consideration 
in spite of the deviation between, in particular, the C values computed in the current study 
and those resulting from the bulk-flow analysis. Otherwise, it is worth noting that the trends 
of the other rotordynamic coefficients are similar in both sets of results and that the direct 
stiffness coefficients, in particular, are in good agreement. Importance of the latter observa- 
tion stems from the fact that the two computational modes are vastly different in the aspects 
of strategy and numerical details, 

4. CONCLUDING REMARKS 

The analysis outlined in this paper provides a detailed look at the rotor/fluid interaction 
in seals as the rotor operation experiences a conical type of whirl. This unstable operation 
can be initiated by a rotor misalignment, as is frequently the case in turbomachinery 
applications, or by rotor self excitation. Regardless of the instability origin, it is crucial to 
identify the sources and magnitudes of the destabilizing moments for the purpose of 
suppressing this type of fluid-induced vibration. 

Applicability of the current finite element-based model is illustrated through a sample 
case involving a straight annular seal to which a simpler bulk-flow model was also applied. 
Comparison of the numerical results in this case led to favorable conclusions, especially in 



Moment coefficients of incompressible-flow seals 


167 


the category of stiffness coefficients. Yet the authors realize that the power of the new model 
would be most appreciated as the model is applied to more complex secondary passages in 
turbomachines, such as that in Fig. 3, where the flow field in the rotor-to-housing gap is 
significantly more complex by comparison. Extension of the model to these and other 
traditionally difficult configurations is currently under consideration. 

Acknowledgement — This study was funded by NASA-Marshal! Space Flight Center (Huntsville, Alabama), 
Contract no. NAS8-37821, NASA monitor James Cannon. Funding was also provided by Houston Advanced 
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A comprehensive approach for computing the dynamic coefficients of an annular 
seal is presented. The coefficients are partly those associated with a uniform lateral 
eccentricity mode of the rotor (known as the cylindrical whirl mode) and with an 
angular eccentricity (which gives rise to a conical whirl type). The rotor excitation 
effects in both cases are treated as interrelated by recognizing the fluid-exerted 
moments resulting from the lateral eccentricity and the net fluid force resulting from 
the angular eccentricity. In all cases, the rotor is assumed to undergo a whirling 
motion around the housing centerline . The computational procedure is a finite- 
element perturbation model in which the zeroth-order undisplaced-rotor flow so- 
lution in the clearance gap is obtained through a Petrov-Galerkin approach. Next , 
the rotor translational and angular eccentricities, considered to be infinitesimally 
small , are perceived to cause virtual distortions of varied magnitudes in the finite 
element assembly which occupies the clearance gap. Perturbations in the flow vari- 
ables including , in particular , the rotor surface pressure, are then obtained by 
expanding the finite-element equations in terms of the rotor eccentricity components . 
The fluid-exerted forces and moments are in this case computed by integration over 
the rotor surface, and the full matrix of rotordynamic coefficients, in the end, 
obtained. The computational model is verified against a bulk-flow model for a 
sample case involving a straight annular seal. Choice of this sample model for 
validation was made on the basis that no other existing model has yet been expanded 
to account for the mutual interaction between the cylindrical and conical rotor whirl, 
which is under focus in this study. 


Introduction 

Annual seals of the type used in turbomachines are known 
to have a significant impact on the rotordynamic behavior of 
the entire engine. Much of the research effort in the turbo- 
machinery community has therefore been focused on creating 
computational tools to predict the stiffness, damping and in- 
ertia coefficients of these seals (e.g., Childs, 1982a; Childs, 
1982b; Childs and Kim, 1985; Nelson, 1985; Dietzen and Nord- 
mann, 1988; and Tam et al., 1988). Simultaneously, several 
investigators (e.g., Childs and Kim, 1985; Childs and Kim, 
1988 as well as Kanemori and Iwatsubo, 1989) have successfully 
measured these coefficients over a range of operating condi- 
tions. 

Perhaps the most widely used model of seal rotordynamics 
today is the Childs* bulk-flow model for its simplicity and 
flexibility. This model was progressively modified to include 
the effects of fluid prerotation (Childs, 1982a), surface rough- 
ness (Childs and Kim, 1985) and seal taper (Nelson, 1985). The 
model was also adapted to handle the case of conical rotor 
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whirl (Childs, 1982b) which results from a tilting (angular) 
eccentricity of the rotor axis. The model is therefore exclusively 
capable of predicting the mutual interaction effects of cylin- 
drical and conical whirl, namely the fluid-exerted moments 
associated with a uniform lateral eccentricity of the rotor axis, 
and the fluid forces associated with a tilting eccentricity of the 
axis. This is the problem under focus in the current study 
(Fig.l), except that the present analysis is much more detailed 
and applies to virtually any leakage passage configuration that 
may give rise to a highly complex flow behavior (Fig. 2). 
Shortcomings of the bulk-flow model, by comparison, emerge 
from a parabolized leakage flow assumption, which would be 
inapplicable should flow separation and recirculation occur, 
and unsimulated fluid stresses away from the rotor and housing 
surfaces. 

Recently, successful attempts to utilize sophisticated com- 
putational fluid dynamics techniques to the seal rotordynamics 
problem have been published. These are best represented by 
the study of Dietzen et al. (1987), in which a traditional per- 
turbation approach to the problem was undertaken, where the 
zeroth-order flow field was computed by solving the time- 
averaged Navier-Stokes equations, together with a two-equa- 
tion turbulence closure, in the clearance gap. However, this 
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Fig. 1 Translational and angular rotor accantricities in annular seals 


computational model, as well as other similar models, have 
not yet been applied to any rotor eccentricity mechanisms other 
than the simple uniform displacement, under which the rotor 
whirl is strictly cylindrical. Other existing detailed models in- 
clude those by Dietzen et al. (1988) and Tam et al. (1988) where 
the flow is instead analyzed in the physically distorted rotor- 
to-housing passage. These are highly accurate models which 
further advocate the use of rigorous flow computational tech- 
niques in predicting the fluid/rotor interaction effects. How- 
ever, such models would be hampered the most by geometrical 
complications should the rotor tilting eccentricity be super- 
imposed on the uniform lateral displacement. 

A common feature in all existing perturbation models of the 
problem is that the perturbation equations (which emerge from 
the flow-governing equations in their differential form) contain 
the clearance gap width. This would imply that the compu- 
tational domain is necessarily bound by the rotor surface, being 
the inner boundary, and the housing surface as the outer 


boundary, and that this condition will have to be valid at all 
axial locations. The deficiency here is that a general leakage 
passage, such as that in Fig. 2, which is naturally connected 
to the primary flow passage at both ends to facilitate primary/ 
secondary flow interaction in real turbomachines (Baskharone 
and Hensel, 1989) would not constitute an acceptable com- 
putational domain for any of the existing perturbation anal- 
yses. The current perturbation model, on the other hand, utilize 
what may be referred to as the ‘‘virtually distorted finite ele- 
ment” technique, whereby the perturbation equations are de- 
rived from the discrete finite-element equations associated with 
the finite elements occupying the flow domain. The latter in- 
cludes, but is not restricted to, the annular passage within which 
the rotor whirl occurs. 

The assumptions made in the present analysis give rise to 
specific applicability limitations which are typical in pertur- 
bation models of the fluid/structure interaction phenomenon. 
For instance, the analysis is inherently inapplicable to large 
eccentricities of the rotor axis (Fig. 3). Moreover it is assumed 
that the clearance-gap flow Field, relative to a whirling frame 
of reference that is attached to the rotor (Fig. 4), is steady at 
all times. This, in effect, rules out seal configurations where 
the rotor and/or the housing are longitudinally or spirally 
grooved. It also precludes situations where the relative flow 
field in the clearance gap is influenced by time-dependent 
boundary conditions at the flow inlet and exit stations. Never- 
theless, the current model embraces the seal category where 
the rotor or the housing surface is circumferentially grooved 
(e.g., labyrinth seals) and tolerates a pre-existing rotor eccen- 
tricity, which is normally a result of misalignment. 

Analysis 

The rotor virtual eccentricity under focus is shown in Fig. 
1. This consists of a uniform lateral eccentricity “e” on which 
an angular displacement “0” is superimposed. These infini- 
tesimally small deviations, from the centered-rotor position, 
are simultaneously accompanied by a rotor whirl at a finite 
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Fig. 2 Generalization of the fluid/rotor Interaction problem: Definition 
of the computational domain and example of the flow behavior In the 
shroudtohouslng passage of a centrifugal pump 



a) CYLINDERICAL WHIRL b) CONICAL WHIRL 

Fig. 3 Cylinderlcal and conical whirl components of the rotor eccentric 
motion 


frequency fl (Fig. l).The fluid reaction in this case is in the 
form of forces and moments in the “jr” and **y” directions, 
which can be related to the displacement components of the 
whirling rotor (Childs, 1982b) as follows: 
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fluid/rotor system, respectively. A single subscript of these 
coefficients implies force-displacement or moment-tilt types 
of relationship, while double subscripts tie fluid-exerted mo- 
ments to rotor displacements and forces to rotor tilting angles. 
Note that a rotor eccentricity, be that linear or angular, in the 
x-direction (for instance) will generally produce a restoring (or 
aggravating) fluid reaction in the ^-direction. Also note that 
straight annular seals are incapable of displaying any cross- 
coupled inertia coefficients, i.e. 

m = m e = m (e = m 0( - 0 (2) 

The procedure to compute the rotordynamic coefficients in 
the matrix equation (1) is initiated by expressing the local radius 
(5) of the rotor axis orbit (Fig. 1) as follows: 

6 = e + e e (3) 

where: 
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( 1 ) 


where the symbols (K,k), (C,c), and (M,m) are the direct and 
cross-coupled stiffness, damping and inertia coefficients of the 


te=0z (4) 

with e e being the lateral eccentricity due to the rotor tilting 
motion. Motion of the rotor axis in Fig. 1 can be decoupled 
into the cylindrical and conical whirl modes depicted in Fig. 
3. Referred to the frame of reference in Fig. 3, the lateral 
displacement components “X" and “7” of the tilting motion 
center can be written for the whirling axis at any given time 
‘7” as follows: 
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Fig. 4 Lateral motion of the coordinate axes before and after the rotor 
eccentricity 


jT=ecos(0/)t y=€sin(Q/) (5) 

Similarly, the angular displacement can be expressed in terms 
of the whirl frequency (Fig. 2(b) as follows: 

e x = - cos(n/), 6 y = - sin(Qr) (6) 

z z 

With no lack of generality, consider the special location of 
the displaced rotor axis where X=0 and 6 y = 0 (as shown in 
Fig. 3). Using expressions (5) and (6) above, the linear and 
angular magnitudes of the rotor axis position, velocity and 
acceleration of the rotor axis at the axial location can then 
be written as follows: 


*=o, y=e 

(7) 

X=Qe, y= 0 

(8) 

su 
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II 
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ll 
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e x =o, 

(ii) 

e x = -tfe, e,=o 

(12) 


Substituting (7) through (12) in (1), and taking the limit as the 
rotor eccentricity components tend to zero, the following equa- 
tions are achieved: 


-k-QC 

3e 


3F i 
de 
dF> 

36 

36 

3M } 

3e 


= -K+Qc + 

— — GC^ + 0 2 

- K ( e + Oc^ 

- _ - GC* 


= K$ t - Gc* e - G 2 M 9( 


(13) 

(14) 

(15) 

(16) 

(17) 

(18) 


in which the cross-coupled inertia coefficients were set equal 
to zero (equation (2)). 

Determination of the rotordynamic coefficients in the dif- 
ferent matrices of equation (1) now reduces to computing the 
derivatives on the left-hand side of equations (13) through (20) 
at a minimum of three whirl frequencies. Interpolation of these 
derivatives as parabolic expressions of 0, using curve fitting 
techniques, leads to these coefficients by simply equating the 
different terms in these expressions to the right-hand sides of 
equations (13) through (20). With this in mind, the remainder 
of this section is therefore centered around computing the 
derivatives appearing on the left-handed side of equations (13) 
through (20) as the ultimate goal of the finite-element pertur- 
bation model. This model, as seen next, is initiated by solving 
the zeroth-order flow field in the undisplaced-rotor clearance 
gap. 

Centered-Rotor Flow Field. This zeroth-order flow problem 
is defined in a frame of reference that is attached to the rotor 
and is rotating (Fig. 4(a)) at the rate of G. Although the axes 
rotation here may appear as serving no purpose at this stage, 
it ensures compatibility of the unperturbed and perturbed flow 
fields, with the rotor eccentricity effects being the only dif- 
ference. As seen in Fig. 4(6), definition of the axes in this 
manner eliminates the time dependency of the perturbed flow 
field once the rotor undergoes the eccentric operation mode. 
In this case, the coordinate axes become of the rotating-trans- 
lating (or simply whirling) type (Fig. 4(£)), with the magnitude 
of Q being set equal to the rotor whirl frequency. However, 
the term G is viewed as an arbitrary rotational frequency of 
the coordinate axes in formulating the undisplaced-rotor How 
problem (which is under consideration at this point). 

Denoting the relative velocity components in the rotating 
frame of reference (of which the origin is fixed, Fig. 4(a)) by 
w, u, and w y the momentum and mass conservation equation 
in the undistorted clearance gap are: 
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(24) 


where: 

u,u,w are the relative velocity components. 

p is the static pressure. 

Q is the rotational frequency of the coordinate 

axes (to be set equal to the whirl frequency 
as the rotor enters its disturbed operation 
mode) 

v h v t> v e are the molecular, eddy and effective com- 
ponents of kinematic viscosity coefficients, 
respectively, where 

V e =Pi+V t 


with the symbol ( " ) designating values that are carried over 
from a previous iteration or an initial guess for the purpose 
of successively linearizing the momentum equations. As seen, 
the whirl frequency is part of these equations, which is a 
result of including the centrifugal and Coriolis acceleration 
effects produced by the axes rotational motion. 

Simulation of the flow turbulence is based on the algebraic 
eddy- viscosity turbulence model by Baldwin and Lomax (1978). 
However, an adjustment in this vorticity-based model was 
implemented whereby the relative, as opposed to the absolute, 
vorticity was used to calculate the mixing length in the fluid 
layer that is adjacent to the rotor surface. Equally important 
in implementing the model was the analysis of the near-wall 
zone which was based on the approach by Benim and Zinser 
(1985). 


Boundary Conditions 

These include known profiles of the inlet velocity compo- 
nents and zero streamwise diffusion of the exit velocity vector. 
The latter boundary condition is achieved by equating to zero 
the second streamwise derivatives of all velocity components. 
This condition is consistent with the foundation of the bound- 
ary layer theory, provided that the boundary layer remains 
attached to the rotor and housing surfaces as the exit station 
is approached. Note that the traditional zero velocity-gradient 
exit condition would be less accurate by comparison, for it 
implies a fully developed flow that may or may not prevail at 
the exit station depending on the seal length. 

The last category of boundary conditions concerns the rotor 
and housing surfaces, and is largely influenced by the rotation 
of the coordinates axes. Referring to Fig. 4 (o), the housing 
surface, to an observer in the rotating frame of reference, will 
no longer appear stationary, but will rather possess a relative 
velocity component Qb as indicated in the Figure, where 
is the housing radius. Furthermore, the rotor surface will ap- 
pear to the same observer as rotating at an angular speed that 
is less than the rotor operating speed a> by the amount f L The 
rotor surface velocity in this case is (co - fi)<7, where “a” is the 
rotor radius. 

Finite-Element Formulation. Numerical solution of the cen- 
tered-rotor flow field is based on the Petrov-Galerkin weighted- 
residual principle in which the weight functions are so chosen 
to incorporate the upwinding concept. An effective tool in 
alleviating spurious pressure modes, this concept has been ad- 
vocated by many analysts in handling highly convective flows, 
such as the present. The concept is equally opposed by others 
due to the false “numerical” diffusion it imparts to the finite 
element equations. The desire to minimize this false diffusion 
in the current flow model led to selecting a mixture of the 


weight functions devised by Heinrich and Zienkiewicz (1977) 
and those used in the conventional Galerkin method in such 
a way that upwinding is applied only to the nonlinear con- 
vection terms on the left-hand side of the momentum equations 
(21) through (23). Note that excessively narrow and long clear- 
ance gaps in typical annular seals would naturally lead to finite 
elements with unavoidably large aspect ratios, especially near 
the walls, due to the special refinement of elements that is 
required across the boundary layer. Nevertheless, utilization 
of the preceding upwinding technique gave rise to a sufficiently 
accurate flow field for elemental aspect ratios of as large as 
196 in the immediate vicinity of the rotor and housing surfaces. 

The finite-element model is constructed with the twenty-node 
quadratic element in Fig. 5 as the discretization unit. The 
procedure is initiated by defining the element shape and the 
cartesian-to-local system mapping relations. Within a typical 
element, let the spatial coordinate conversion be defined, in 
terms of the nodal values, as follows: 

20 20 20 

X = jytiKi , y = J^N.yi , z = (25 ) 

/.] 1*1 /-l 

where TV/ are quadratic “shape” functions (Zienkiewicz, 1971) 
of the local coordinates f, ij and { (Fig. 5). Conversion of the 
spatial derivatives is, in this case, defined as follows: 



where the components of the matrix [T] are as follows: 



with I J 1 being the Jacobian of transformation. 

Next, the flow variables are interpolated within the typical 
element (Fig. 5) in a similar fashion. Guided by the Ladysh- 
enskaya-Babuska-Brezzi compatibility requirements (Carey and 
Oden, 1986), as applied to the current problem, the velocity 
components and pressure are expressed as follows: 


v 
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Fig. 5 Quadrilateral curve-sided Isoparametric finite element for ana- 
lyzing the rotor-to-housing flow field 


20 20 20 8 

u = y/V.u, , v = ]>>, v t , w = y]N, W, , p= ^TJi MkPk ( 28 ) 

#- 1 r's] i-1 

where the interpolants “A/*” appearing in the pressure expres- 
sion are the linear shape functions associated with the finite 
element corner nodes (Zienkiewicz, 1971). 

The error functions produced by the unperturbed flow gov- 
erning equations, as a result of the interpolation expressions 
above, are then made orthogonal to a special set of weight 
functions W x throughout the finite element. In constructing 
these functions, the so-called error consistency criterion of 
Hood and Taylor (1974) is implemented, whereby the element 
linear shape functions, A//, are used in conjunction with the 
continuity equation. The weight functions used with the mo- 
mentum equations, on the other hand, are of the upwind form 
described by Heinrich and Zienkiewicz (1977) when applied to 
the convection terms, and are identical to the element quadratic 
shape functions ‘W/’ for all other terms. 

The orthogonality conditions, stated above, constitute the 
elemental set of equations, which can be written as follows: 

[j,] ( u ) + [a 2 ] M + [a 3 ] [ w ) + [a 4 ] [p ) = [ b, ) (29) 

fol ( * ) + fo] ( v ) + lail ( w ) + [a 8 ] [p } = { b 2 ) (30) 

[a 9 ][u] +[aiaHv) +[flu]M + = [b 2 ] (31) 

[013]|H] + [*m]M + faisl M =0 (32) 

where the vectors {«), [ t? } , (w), and (/?) contain the nodal 
values of the velocity components and static pressure, asso- 
ciated with the finite element. The coefficient matrices [aj] 
through [tfj 5 ] and the right-hand side vectors (/?i) through 
[6 3 ) are all defined in Appendix I. Equations (17) through 
(20) can be rewritten in the following compact form: 

[*](*} = [£) ( 33 ) 

where the vector (<£} contains the nodal values of velocity 
components and static pressure associated with the typical 
finite element. The global form of the elemental set of equa- 
tions (33) is obtained by assembling all finite element contri- 
butions, and introducing the various boundary conditions, at 
which point the equations are expressible in the following form: 

H](*) = {2?) (34) 

where the global vector ($) contains all nodal values of the 
velocity components and pressure prior to the rotor excitation. 
These unperturbed flow variables constitute an axisymmetric 
flow field for a perfectly centered rotor, which is the special 
case under consideration in the current study. However, should 
a rotor misalignment (including a fixed tilt) be the case, the 
unperturbed flow field described by equation (34) would gen- 
erally be three-dimensional and should be modeled as such. 

Perturbation Model. Referring to equations (13) through 
(20), determination of the rotordynamic coefficients reduces 


to the fundamental problem of computing the rate at which 
the fluid-exerted forces and moments in the clearance gap are 
developed with respect to a compound (translational and tilt- 
ing) eccentricity of the rotor axis. This problem is analyzed 
next. 

The current computational model is centered around the 
manner in which the finite-element equations (34) are altered 
as a result of the rotor virtual eccentricity depicted in Fig. 1. 
Since the rotordynamic coefficient are (by reference to equa- 
tions (13) through (20)) dependent on partial derivatives of the 
rotor displacement components “e” and **0", it is legitimate 
to treat the fluid response to each component separately. Re- 
ferring to Fig. 3 in this case, let "e” symbolize the rotor axis 
eccentricity, where: 

e for translational motion ^5) 

6z for tilting motion 

with z being the distance from the center of the tilting motion. 
Under such displacement, each finite element will, as a result 
of virtual distortion (Fig. 6), yield a perturbed set of equations. 
Let the assembled form of these equations, by reference to 
(34), be as follows: 

([A] + e[A))(i*) + (6$!) = ([5) +e{ B}) (36) 

Subtracting (34) from (36), the following expression for 
is obtained: 

(5$j=^]- , ((S)-Ml(4>)) (37) 

or, 

= (38) 

be f-o\e/ 

The general equation (38) implies that two new arrays, namely 
the matrix [A] and the vector { B ) , are needed to compute the 
rate at which the flow variables (namely the velocity compo- 
nents and pressure) vary with respect to the rotor eccentricity 
"e.” These arrays are primarily a result of finite elements 
distortion effects that are associated with the rotor displaced 
position (Fig. 6). The arrays [A] and (£) are derived next as 
part of the procedure to compute the rotordynamic coeffi- 
cients. To do this, a typical distorted element will first be 
analyzed, and the elemental contributions [a] and ( b 1 derived. 
The global arrays [A] and { B j can subsequently be obtained 
by assembling [a] and ( b ) among all finite elements. 

Analysis of a Distorted Finite Element . Consider the dis- 
tortion experienced by the typical finite element shown in Fig. 
6. As seen, the distortion is caused by virtual displacements 
of the element nodes as a result of the rotor eccentricity. Re- 
calling that the axes are attached to the whirling rotor, it is 
those nodes on the housing that will be displaced the most, 
while the rotor surface nodes remain stationary. 

With “e” signifying the linear or angular eccentricity 
(expression 35), let the new coordinates of a typical node 
(Fig. 6) be (X iy Y it Zj), where: 

Xi-Xi, Y, = y t + X , e and Z t = z, (39) 

where in general: 

X, = M Xi,y it Zi) (40) 

Dependency of X, on the in-plane position of the node **i r 
(i.e., Xi and y,) is clear by reference to Fig. 3. However, it is 
only the purely conical whirl in this figure where X, is addi- 
tionally dependent on “Z/ For instance, X, = -z, for nodes 
on the housing surface under this type of whirl (note the po- 
sitive direction of the y-axis in Fig. 6). 

Let the new local-to-cartesian transformation Jacobian 1 Ji I 
be related to the undistorted-element Jacobian I J I as follows: 

I Jj I = I JI +el J I (41) 

The determinant 1 J I can be derived using the coordinate re- 


Joumal of Tribology 


JULY 1991, Vol. 113/475 



-6F , 


Fig. 6 Deformation of the finite elements In the roto-to-houslng gap 
as a result of the rotor eccentricity 


(21) through (24) and reapplying the weighted-residual pro- 
cedure, the distorted element equations can be written (upon 
separation of terms containing e) as follows: 

(W+e[fl)XI«l + W))*(l*}+e|6)> (45) 

where the matrix [a] and vector [ b ) result from the distortion 
in the element shape. The vector ( 50 ) in this case contains the 
differential changes in the nodal values of velocity components 
and pressure that are created by the virtual eccentricity. 

Noting that the matrix [a] is similar in construction to [a] 
in equation (33), the former can consistently be viewed as 
composed of submatrices [flj] through [o I5 ], which are com- 
patible with those in Appendix I for the undistorted finite 
element. Jhese new matrices are defined in Appendix II. The 
vector (b) is substructured in a similar manner, and is also 
defined in Appendix II. Finally, the global arrays [A] and 
[B] appearing in equation (38) are attained by assembling 
[a] and { b ) among all finite elements. 


lationships (expression 39) and ignoring high-order terms of 
“e” as follows: 

k ^ Xi h ‘ h z< 

yiM, y>9N, ys dN, 

kn x ‘ kn Zi 

A new spatial derivative conversion relationship, that is equiv- 
alent to equation (26) can now be derived as follows: 


* ) = ([71+e[P1) i a 


(43) 


where the matrix [T] is as described in equation (27) while the 
matrix [P] is defined as follows: 

p = _L ( y >dN ‘ x V— 7 ~ y ,dNi x V 1 — -ill t 

" I J I \,T{ <>v ‘h ^ h K ‘h 7 I J 1 2 " 

P -J-fyMi YdN, fa*; \ m 

12 iji V"d{ kn k^ Xi kn 7 ui 2 12 

P = _L /f x f^ 7 _ f m x f m , V ill T 

13 ui \k d t x, ,r|at) ' k d 7 x 'r ! ^ / ui 2 7,3 


„ ui 

P21 TjF 721 

„ i j i „ 

P22= "TjP Tu 

n I J I „ 

Pa= ~ iTP 723 
20 


(44) 


p = _L (y™i x - x\-UL r 

Pjl IJI dr, k ^ X ' k K k dl > 7 UI 2 731 

P = _L xV— X-Y 1 — xrV'— T 

32 ui \kn kn 1 k ^ k n 7 ui 2 32 

p = -L(v™i x y™i x-fe x^-ilir 

33 ui v& 3 f k ‘ k ** h ^ 7 ui 2 33 

Substituting expressions (41) and (43) into the flow equations 


Calculation of the Rotordynamic Coefficients . Equation (38) 
can be rewritten in the following detailed form: 

j-w-aii-wi*,) ( 46 ) 

where ( U} t { K), { W), and (P) are global vectors containing 
the nodal values of the velocity components and pressure, 
respectively, with the overbars signifying quantities that are 
associated with the perturbation in the rotor-to-housing flow 
field. Of all the pressure nodal values in the vector (P) above, 
define a subvector ((P) as composed of the pressure values at 
the rotor surface nodes, i.e., 


d_ 

de m de 



{<P) = (A, i=W) C{P) (47) 

where N s is the total number of comer nodes existing on the 
CdfP) /(dp\ 1 

rotor surface. Note that j— j = M— ), /=1,AM is but a 


subvector of the already computed global vector — in equa- 




tion (46) and is, therefore, known at this computational step. 

Next, consider a finite element face “ s M that exists on the 
rotor surface (Fig. 7). The pressure derivative over this 
surface element can now be interpolated in terms of the known 
pressure derivatives at the four corner nodes as follows: 

(48) 

Once integrated over the proper differential segments of the 
rotor surface, the function p t [ s) (above) yields a differential 
force which, in turn, can be resolved in the “jt” and 
directions. Furthermore, taking the moments of these forces 
around the center of tilting motion in Fig. 3(b) and summing 
the result over all finite elements sharing faces with the rotor 
surface, the rate (with respect to “e M ) at which the fluid re- 
action moments are exerted can be attained. The different 
derivatives in equations (13) through (20) are then obtained 
by appropriately replacing e with e or dz (equation (35)) as 
follows: 

dF N$ f +1 f +1 

-^ = JC J - , J - 1 (49) 

dF ^ |* + 1 + 1 

}_,nJv.OP„(v.S)G(v.WvdS (50) 

dAf ^ f + * f * * 

]_| znjXi),()p, f (r) f ()G(ri t ()dridZ (51) 

j * i 
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Fig. 7 integration of the pressure forces on the rotor surface 


d A/ * /* + i /• + 1 

-^=2 j-i j-i zn/n^)p, t (v,i)G(ri,i)dvdi (52) 
dF ^ 5 f + * [* + * 

^=E]_, J . [ Z n A*l,Z)P>i)(.v,£)G(y,£)dr)dt (53) 

dFy ^.PT' 

-^= L J- ! J - 1 znAvMPMMvMdvdi (54) 

f « 1 

dM y v^p'p 1 

l/ = LJ-, )- [ ™M,Z)P,e(ri,k)G(r,,i)d V di. (55) 
dM x ^f + 1 f + l 

-^"=L J-, (56) 

where n, and n, are the components of the local unit vector 
that is normal to the rotor surface. Also, the parameter G(i j,{) 
is a Jacobian-like function for cartesian-to-local area trans- 
formation, and was previously derived by Baskharone (1979) 
as follows: 

g( t )=\ (—— V 

LU;3f dtdv) + W d n d v dj 

dx dy dxdyVli 

J 

where the derivatives on the right-hand side are evaluated using 
the interpolation expressions (25) for a “f” value of 1.0 (Fig. 
7). As seen, the first four derivatives (equations (49) through 
(52)) are the outcome of a cylindrical whirl problem (Fig. 3(a)) 
where e = e, while the remaining derivatives are associated with 
the conical whirl problem, in which case e = Bz. (Fig. 3(b)). 

At this point, the requirements for computing the rotor- 
dynamic coefficients in equation (1) are, by reference to equa- 
tions (13) through (20), complete. The process, as described 
earlier, consists of computing the force and moment derivatives 
in equations (49) through (56) at three or more whirl frequen- 
cies. Substitution in equations (13) through (20) finally leads 
to the different rotordynamic coefficients by inspection. 

Results and Discussion 

A hydraulic annular seal with a clearance/length ratio of 
0.0085 was selected for the purpose of verifying the current 
computational model. The force coefficients of this seal were 
determined by Dietzcn and Nordmann (1987), and later verified 
by Hensel (1990). This seal has a rotor radius of 23.5 mm and 
a Reynolds number (based on the clearance and inlet through- 
flow velocity) of 4700. Operating speeds of 2000, 4000, and 
6000 rpm were considered in this case, together with preswirl/ 
inlet velocity ratios of 0.04, 0.10, and 0.17, respectively. The 
zeroth-order flow solution corresponding to each set of op- 


erating conditions required an average of 3.5 hours of CPU 
time on the VAX 8800 mainframe. Execution of the pertu- 
bation analysis, on the other hand, consumed nearly 1 2 minutes 
of CPU time on the NEC SX2 Supercomputer for each whirl 
frequency of the rotor axis. 

An early version of the current analysis was tested by Hensel 
(1990) using the same annular seal described above, and the 
results compared to experimental and analytical data. Missing 
in these results, however, were the rotordynamic coefficients 
linking the effects of translational and tilting rotor eccentricies 
to one another since the computational model, then, did not 
account for such effects. Presentation and assessment of these 
coefficients are therefore the focus of this section. 

Shown in Fig. 8 is a comparison between the computed values 
of the moment coefficients. K Sty k 6t , C 9f , and c Bi associated 
with the rotor cylindrical whirl, and those obtained through 
Child’s bulk flow model. Of these, the cross-coupled stiffness 
and direct damping coefficients k 9t and C d( represent the mo- 
ments around the y-axis (Fig. 2) which arise from the rotor 
cylindrical whirl. These moments suppress (or aggravate) the 
rotor whirling motion and are therefore of particular interest. 
As seen, the trends of these two coefficients with the rotor 
spinning speed are basically similar in both sets of results. 
However, the bulk-flow model seems to produce magnitudes 
that are significantly smaller than those produced by the cur- 
rent analysis. 

Figure 9 is intended to support, in a general sense, the impres- 
sion that the bulk-flow model tends to underpredict the coef- 
ficients k 9< and C* discussed above. The comparison in this 
figure is between the results of applying this model (Childs, 
1982b) and the experimental measurement by Kanemori and 
Iwatsubo (1989) for a hydraulic annular seal with a clearance/ 
length ratio of 0.002 and a length/inner diameter ratio of 3.0. 
The rotordynamic coefficients in this figure were obtained by 
interpolation at a fixed overall pressure drop of 100 k. Pa. 
Examination of Fig. 9 reveals that this model is constantly and 
significantly underpredicting both coefficients. It should be 
pointed out here that conversion of the experimental data by 
Kanemori and Iwatsubo (1989) into the form shown in Fig. 9 
may have involved some interpolation errors, and that the 
differences (in dimensions and operating conditions) between 
the tested seal and the current, indeed present another source 
of uncertainty. Nevertheless, the differences between our com- 
puted magnitudes of k e< and C 9e (Fig. 8) and those produced 
by the simple bulk-flow model, seem to be consistent with the 
experimental findings. 

The rotordynamic coefficients K (9f k^ y C e9l and c< 6 in Fig. 
10 represent the fluid-exerted forces on the rotor as a result 
of conical whirl. Again, the cross-coupled stiffness and direct 
damping coefficients ( k # and C rf ) give rise to aggravating and 
restoring forces, respectively, in the event where k t9 is negative 
which is the case here. Aside from the moderate agreement 
between the computed values of K& and c (6 and those obtained 
through the bulk-flow model, it is seen that the respective 
values of k& and C ( $ are far apart. In fact, our results indicate 
that the rotor is comparatively much more stable, as far as 
this conical-to-cylinderical whirl conversion is concerned. Un- 
fortunately, there is no experimental or analytical data covering 
this aspect of rotor eccentric motion that can be used to favor 
a set of numerical results over the other, as the case was in 
discussing the counterparts of these coefficients earlier in this 
section. 


Concluding Remarks 

The perturbation model presented in this paper is a quite 
general tool that embraces the most general mode of rotor 
eccentric motion and lends itself to a challenging category of 
fluid/rotor interaction problems, yet to be explored. This is 
where the seal is part of a secondary (or leakage) passage which. 
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Fig. 8 Moment coefficients due to the rotor translational eccentricity 
and comparison with Childs' bulk-flow model results 
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imental data of Kanemori and Iwatsubo (1989) 
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in turn, is connected to the primary flow passage at both ends 
in a manner that permits the mutual interaction between the 
two flow fields (Fig. 2). The authors are not aware of any 
existing perturbation analysis that is conceptually applicable 
in this case for the simple fact that existing analyses address 
the problem where the gap between the whirling body and the 
housing extends between the flow inlet and exit stations, with 
the gap width being part of the flow-governing equations. This 
restriction is totally alleviated in the current model since the 
virtually-distorted finite element concept, on which the com- 
putational model is based, would in this case be utilized ex- 
clusively in subdomains which interface with the whirling rotor 
surface. The difficulty, however, in investigating such a prob- 
lem stems from the need to use a large size finite-element model 
(Fig. 2) as the computational domain is unavoidably complex 
from a geometrical standpoint. Efforts are currently made by 
the authors to minimize the computational resources needed 
in this case towards a sufficiently accurate numerical solution. 

Current analysis of the coupled cylindrical/conical whirl 
effects in annular seals is preceded only by simplified models 
which do not fully account for the flow details in the annular 
passage. The analysis yields the direct and cross-coupled stiff- 
ness and damping coefficients linking each whirling motion to 
the other. Comparison of the results with those obtained 
through Childs’ bulk-flow model revealed that the latter model 
generally underpredicts the rotordynamic coefficients under 
focus in this study. Experimental evidence, however, suggests 
that the current results are more accurate by comparison. 
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APPENDIX I 

Coefficient Matrices of an Undistorted Finite Element 


Define the differential operators F, G, and II such that: 


d d d 

F = 7 -„- + r,:- + r 13 - 

d d d 

G =7' 21 -+r 22 - + r 23 - 

d d d 

H = 7 ', 1 - + r 32 - + r 33 - 


where the coefficients T xu T l2 ,...e tc. are as defined by the set 
of expressions (27). The different coefficient arrays in the 
equations associated with a typical finite element (equations 
(29) through (32)) can now be expressed in terms of these 
operators as follows: 

J,,, I P e [FM)F (Nj) + GmG(Nj) + H(N,)H(N y )] 


- W'[G(f,)G(N / ) + H(P ; )H(N,) + 2F(P,)F(N,)] 
+ H'lwFfA/y) + vG(Nj) + wH(Nj)) \dV 


02),, = NJGftJF(N/) - 2QNj]dV 

a, ),j= J*„- W,H(v,)¥(Nj)dV 
*«),.*= N?{M k )dV 

«5 hi- J *.» - N/[F(?,)G(Ny) + 2QNj]dV 

««),.;= J^r) ( iv[F(A mNj) + G(W;) + H(W;)1 


- Wi\FO,W(Hj) + H(J,)H(N,) + 2G(P,)G(N,)] 

+ fV,[u¥(Nj) + CG(Nj) + vvH(jV,)] }dV 

«7 {^„-N,H(P,)G(N,)rfK 

<*«)<,* = ^G(M t )dV 

a,)u= 

«io),.y= j^,-/V,G(P,)H(N>fK 

an)u = J*,, ( W)F(N ; ) + GW)G(AO) + H(N,)H(/V,)] 

- W'lFfrJFCNy) + G(i>,)G(N,) + 2H(P,)H(/V,)] 

+ WiluTWj) + "GW) + kM(N,)] \dV 
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a,2)6*=J^WA/*)rfK 

«13)*.y = [„ () M t F(^ 

«.5)*.y= { k ,„A/*H(^MK 

*.)/= J s ,e,W(n- VH^ZS+Q 2 

* 2 ), = J,,., W<« • v v)dS + n 2 } K(f) N, (L^') <* k 

*3),= J s «>W<n • V w)dS 

where i- 1,2,. ..,20 t j- 1,2,. ..,20, k= 1,2,. ..,8, the symbol f ) 
refers to a value that is known from a previous iteration or 
an initial guess, dS is the differential element of surface area 
(Fig. 5) and n is the local unit vector normal to the boundary. 
Also, the volume differential dV is as follows: 

dV= IJI dtd v dt 

where I J 1 is the Jacobian of transformation from the local to 
the cartesian frame of reference (Fig. 5) for the undistorted 
finite element. 


APPENDIX II 
Perturbed Set of Finite Element Equations 


+ G (N,)3(Nj) + W)G(N t ) + H(N t )K(Nj) + K(N t )H(Nj)] 

- Nff&dW) + 1 ftW + H(*,)K (Nj) + K(P,)H (Nj) + 2G(*,) 
x J (fy) + 23(p t )G(Nj)] + N'luKNj) + C3(Nj) + wK (Nj)] }dV 
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* 

lyit) 

ff io)i,y— a io)ij- 

[ NjlG(v,yK(Nj) + 3(p t )U(Nj)]dV 


Jy(e) 

®n)ij—®fi)u+ 

[ { wmw) + m)T(Nj) + G(N,)3(Nj) 

J yie) 


+ W)G(Nj) + HW)K(N,) + K W)H(N,)] - NffiWWj) 
+ I(WW + G(*,)J(N y ) + J(*,)G(N ; ) + 2H(S>,)K(N > ) 

+ 2K(P,)H(N,) + NilQUNj) + 03 (NJ) + wK(N,)] } dK 

Oi2)/.*=ar2)t*+ [ - NjK(M k )dV 

Jy«)P 

ii2)kj=ah)*j+ f m (Nj)dV 

Jyi*) 

a u)*, y = *«*.;+ [ M^NjyiV 

Jyie) 

a ti ) k j~ai s ) k j+ \ MMN^V 
Jyie) 


Referring to equation (45), the matrix [c] and vector (Z?) 
associated with the distorted element in Fig. 6 is broken up 
into smaller arrays [ a{\ through [a 15 ] and { b \ ) through ( ) . 
These arrays are similar in construction to those defined in 
Appendix I (for the undistorted element), and can be expressed 
as follows: 

a 1 )u = ar\j + [ { vJF(Ni)l(Nj) + l(N t )¥(Nj) 

J y(t) 

+ G(N,)3(Nj) + W,)G(Nj) + H(Wj) + K(W;)] 

- NilG(v,)J(Nj) + J(P,)G(N,) + H(P,)K(N,) + K(P,)H(N,) 

+ 2¥(v,)l(Nj) + 2I(P,)F(N,)] + Niltil(Nj) + CJ (NJ) + *K(N,)] }dV 

( ^[G(?,)I(N,) + J(P,)F(N y )]rfK 

\j=aikj- ( JV,[H(P,)I(N ; ) + K(P r )F(A^)]^F 

Jyi*) 

04 ), * = *;),*+ ( - N\(M k )dV 

a i ),j = a;), j - \ N,[F(P,)J(N y ) + I(P,)G {NffldV 

Jyie) 

= \ ( ?,[F(N,)I(Ar y ) + \(N,)F(Nj) 

J yie) 


m + 1 / 20 \ 

^ 1 J Wdvdt 

m + 1 / 20 \ 

1 j irfUME 






In these expressions, a variable with an asterisk has the same 
form as that in Appendix I, with the exception that the Jacobian 
I J I is now replaced by I j I as defined in equation (41) where 
1 and j vary from 1 to 20, while k varies from 1 to 8. With the 
operators F, G, and H being those defined in Appendix I, the 
new operators I, J, and K are defined as follows: 


d d d 

l=p "K +p »T v +p »Js 
d a d 

i = p 2)-^+ p xy v + p »-^ 

K= p 3 , ?-t +Pn T +p *h 

df dr} 3 £ 


with the matrix [F\ being that defined by the set of expressions 
(44). The volume differential “dV” is the same as defined in 
Appendix I for the undistorted finite element. 
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Flow Field in the Secondary, Seal- 
Containing Passages of Centrifugal 
Pumps 1 

This paper illustrates the impact of seal configuration on the through-flow leakage 
in centrifugal pumps with shrouded impellers. The flow model is based on the Petrov - 
Galerkin finite element method , and the computational domain permits the primary/ 
secondary flow interaction at both ends of the clearance gap . The model is applied 
to a hydraulic pump with two different seal configurations for the purpose of 
comparison . The computed results show a strong dependency of the leakage flow 
percentage and swirl-velocity retention on the overall shape of the shroud-to-housing 
passage including, in particular, the seal geometry. The results are generally con- 
sistent with documented observations and measurements in similar pump stages . 
From a rotordynamic standpoint , the current computational model conceptually 
provides the centered-rotor tt zeroth-order”f1ow field for existing perturbation models 
of fluid /rotor interaction. The flow model is applied to two different secondary 
passage configurations of a centrifugal pump , and the results used in interpreting 
existing rotordynamic data concerning the same passage configurations. 


Introduction 

Leakage flow in the shroud-to-housing gap of centrifugal 
pumps has significant performance and rotor-integrity con- 
^ sequences. First, it is the leakage flow rate, as determined by 
the through-flow velocity component, which is typically a ma- 
jor source of the stage losses. The swirl velocity component, 
on the other hand, is perhaps the single most predominant 
destabilizing contributor to the impeller rotordynamic behav- 
ior. Control of the through-flow velocity in the clearance gap 
is often achieved through utilization of a tight-ciearance seal. 
^ Suppression of the flow swirl, however, requires careful design 
of the leakage passage and/or the use of such devices as the 
so-called swirl “brakes” (e.g., Childs et al., 1991) or straight- 
ening grooves/ribs in the inner housing surface (e.g., Ohashi, 
1988). Unfortunately, an efficient leakage control device, such 
as the labyrinth seal, may itself trigger an instability problem 
of fluid-induced vibration (Childs and Elrod, 1988). In the 
current study, two seal configurations, comprising part of the 
^ shroud-to-housing passage in a centrifugal pump, are analyzed 
for comparison. In both cases, the impeller geometry and the 
pump operating conditions are identical. 

The current computational tool is an expanded version of 
a finite-element model which was previously proposed by the 
authors (Baskharone and Hensel, 1989), where the mere idea 
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of including primary-flow segments in the computational do- 
main definition was introduced. This was a means of avoiding 
the need for what would otherwise be unrealistic boundary 
conditions at the primary/secondary flow interface. The model, 
initially based on a laminar flow' assumption, has since been 
upgraded by recognizing such aspects as turbulence and inertia 
domination of the flow field. 

The outcome of the current study is tightly linked to the 
rotordynamic stability analysis of shrouded pump impellers. 
This is true in the sense that the numerical results can essentially 
be utilized as the “zeroth-order” flow field in existing per- 
turbation models (e.g., Baskharone and Hensel, 1991a) for 
computing the stiffness, damping and inertia coefficients of 
the fluid/shroud interaction forces, as the impeller axis under- 
goes a whirling motion around the housing centerline. Accu- 
racy of these rotordynamic coefficients was reported by 
Baskharone and Hensel (1991a) to be a strong function of the 
centered-impeller flow solution, which is under investigation 
here. 

Analysis 

Figure 1 shows the primary and secondary flow passages in 
a typical pump stage, and the computational domain under 
consideration. The latter includes two primary-flow segments, 
which render the domain to a double/entry, double/departure 
flow region. Also shown in Fig. 1 are the major features of 
the finite element model in which a biquadratic curve-sided 
element (Fig. 2) is used as the discretization unit. 

Definition of the computational domain in the manner shown 
in Fig. 1 is hardly traditional. Inclusion of two primary flow' 
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Fig. 2 Biquadratic nlne-noded finite element of the Lagrangian type 


segments in this domain is aimed at facilitating, to a reasonable 
level, the primary/secondary flow interaction effects at both 
ends of the secondary passage, without resorting to a much 
larger numerical model, had the impeller subdomain been 
added. Existing computational models, by comparison, treats 
the secondary passage as totally isolated from the primary- 
flow passage (e.g., Childs, 1989). As seen in Fig. 1, inclusion 
of the primary-flow passages (labeled a-b-d-c and e-f-h-g) 
clearly alleviates the need to specify what would otherwise be 
grossly simplified boundary conditions at the two primary/ 


secondary flow interaction locations at both ends of the sec- 
ondary-flow passage. 

It might appear, at first, that the current flow' problem is 
solvable using one of the existing commercial flow codes. How- 
ever the very nature of the computational domain, as a multiple 
entry/departure flow’ region, and some of the corresponding 
boundary conditions (discussed later in this section) are too 
nontraditional for such codes. 


Flow -Governing Equations. The momentum and mass 
conservation laws governing the swirling axisymmetric flow' of 
the incompressible fluid in Fig. 1 can be expressed as follows: 

V, ~7^+ V ir-~ = 

dr dz r p dr 

^_j Vr + ^f (I) 

dr dr r dz or 


v 'ir + v - ^+— = v.^v v,) 

dr ' dz r 


[ <3^r »ctf 

r dr ^ r 


V* & 


v dh +v d -±-J-*P 

' dr : dz p dz 


dv t dV z dv, dV r 

+ V.(.,„VK) + — — — (3) 


dr dz r 


(4) 


where: 


V n V d , and V. are the r, 0, and c velocity components 
p * is the static pressure, 
p = is the flow density, 

v t and y C ff = are the eddy and effective kinematic viscosity 
coefficients, respectively. 


Turbulence Model. The turbulence closure in this study is 
that devised by Baldwin and Lomax (1978), together with an 
enhanced version of the near-wall zone treatment proposed by 
Benim and Zinser (1985). First, the effective kinematic viscosity 
i/ efft in Eqs. (1) through (3), is cast in terms of the molecular 
and eddy components, v ( and v t . In calculating the eddy com- 
ponent, p„ the procedure assumes the presence of two, inner 
and outer, layers. In the inner layer, the Prandtl-Van Driest 
formulation yields the following expression: 

p,. , = / 2 lwl (5) 

where the subscript / refers to the inner layer, and the symbol 
u> stands for the local vorticity. The mixing length, /, in expres- 
sion (5) is defined as follows: 
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Fig. 3 Near wall refinement for accurate prediction of the eddy vis- 
cosity 
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v; and 


( 6 ) 


where: 

y = the distance normal to the nearest wall, 

A ’ = the sublayer thickness, 
r w = the w-all shear stress 

The model switches from Van Driest formulation to that of 
the outer region at the smallest value of y for which the inner 
and outer values of the eddy kinematic viscosity are equal. The 
formulation for the outer layer is given by: 


v i. o — K^cpE mixyrmx^KLEB 


( 7 ) 


where: 



with y miy referring to the value of y at which F max occurs. The 
various constants in Baldwin-Lomax model are as follows: 


A* = 26, * = 0.4, * = 0.0168, C f ,= l. 6; and C KLEB = 0.3 

A modification to the Baldwin and Lomax turbulence clo- 
sure was made in the current study to accommodate the angular 
motion of the shroud surface (Fig. 1). Referring to the defi- 
nition of the inner-layer kinematic eddy viscosity, and consid- 
ering the case in which the layer is attached to this rotating 
surface, the vorticity and wall shear stress were both defined 
on the basis of the relative velocity (W) as opposed to the 
absolute velocity (V), where: 

W = V - Qre H (8) 

where: 


= w - 2Q 


(9) 


Calculation of the eddy viscosity in the outer layer (Eq. (7)), 
as pertaining to the rotating surface, was consistently per- 
formed using the relative flow properties. 

Computation of the wall shear stress, t w , is based on the 
near- w all zone treatment proposed by Benim and Zinser ( 1 985). 
The assumption here is that the universal law of the wall at 
any wall location is extendible to an interior computational 
node that is closest to the wall at this location (Fig. 3). Referring 
to this minimum distance from the wall by y min , the following 
expression for the wall shear stress is then obtained: 


^ v lP ^min 
^mm 

Tw ~ r>\ '4 rj/ i,l/2 

1 P rr mm A. min 

: In 

where: 


for 11.6 


for >w,> 11.6 


(10) 


*n,,n = -£h. C D = 0.09, k = 0.4, E= 9.0 

The symbol fF min in expression (10) refers to the interior-node 
velocity relative to the wall for a rotating wall segment, and 
is identical to the absolute velocity otherwise (Fig. 3). Note 
that the outcome of this equation in the case where y^w ^11-6 
is a recursive relationship since the wall shear stress r H now’ 
appears on both sides of the equation. An iterative procedure 
is executed, in this case, to compute r w . 

Enhancement of the accuracy of the preceding turbulence 
closure, including the near-wall model, is achieved with the 
aid of an array of points that is different from the primary 
set of computational nodes in the finite-element model. This 
is in contrast to the model by Benim and Zinser, where the 
flow' properties at the finite-element node closest to the wall 
were used to calculate the wall shear stress. Figure 3 shows an 
enlarged segment of the computational domain near a generally 
rotating wall, in which the primary nodes in the finite-element 
discretization model are identified by hollow circles, while the 
points used in the eddy viscosity computations at the typical 
node “f * are solid circles. The objective here was twofold; to 
estimate the cut-off location between the inner and outer layer 
with sufficient accuracy, and to capture the steep gradients of 
the flow variables near the solid wall. 


Boundary Conditions. Referring to the flow-permeable 
boundary segments in Fig. 1, the boundary conditions over 
these segments are as follows: 

(/) Stage Inlet Station. This is the boundary segment a-b 
in Fig. 1, which is located sufficiently far upstream from the 
impeller. Fully developed flow is assumed at this location, 
giving rise to the following boundary condition: 

dz dz dz ~ 

In addition, the stage-inlet static pressure is specified at the 
node midway between the endwalls on this station. 

(//) Impeller Inlet and Exit Stations. These are labeled c- 
d and e-f in Fig. I. Fixed profiles of the velocity components, 
corresponding to the stage operating conditions, are imposed 
over these boundary segments. Note that the operating con- 
ditions here involve the primary impeller passage, and do not 
include the secondary mass flow rate. 


Q = the rotor spinning speed 
r = the local radius 

e e = the unit vector in the positive tangential direction 
In this case, the relative vorticity (u R ) is expressible in terms 
of the absolute vorticity (u>) as follows: 


(Hi) Stage Exit Station. The flow behavior at this station 
(designated g-h in Fig. 1) is viewed as predominantly confined 
to satisfying the mass and angular momentum conservation 
equations in a global sense. In their derivative forms, these 
can be expressed as follows: 
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These two boundary conditions are linear and are, therefore, 
introduced non-iteratively in the numerical solution process. 
Moreover, a zero normal derivative of V. is imposed over this 
station, and the stage-exit static pressure is specified at the 
computational node midway between the endwalls on this sta- 
tion. 

As for the solid boundary segments in Fig. 1. namely those 
of the housing and shroud as well as the hub surface segment 
b-d, the no-slip boundary condition applies as follows: 

V, = 0, K : = 0 and V e = C 

where Cis equal to (Qr) and zero for rotating and nonrotating 
boundary segments, respectively. 


Finite-Element Formulation. A special version of the Pe- 
trov-Galerkin weighted residual method is used to derive the 
finite-element form of the flow-governing equations. The cur- 
rent approach ensures upwinding of the convection terms in 
the momentum equations while preserving the elliptic nature 
of the diffusion terms. This, for a simple orthogonal grid, 
would be equivalent to backward-differencing the convection 
terms and central-differencing the diffusion terms in the con- 
ventional finite difference analyses of inertia-dominated flows. 
Successful implementation of this strategy, within a finite- 
element context, was achieved by Hughes (1978) for only simple 
(linear and bi-linear) finite-element configurations by modi- 
fying the integration algorithm in the process of deriving the 
element equations. This effectively eliminated the “wiggles” 
in the streamwise pressure variation which are typically as- 
sociated with the conventional Galerkin’s weighted residual 
approach when applied to high Reynolds number flows, such 
as the present. Expansion of essentially the same concept to a 
highly accurate biquadratic element, is employed in the current 
model, by selecting the weight functions on a term-by-term 
basis. 

The characteristic features of the finite element discretization 
model in the current study are shown in Fig. 1. The discreti- 
zation unit, which is a nine-noded curve-sided finite element 
of the Lagrangian type (Zienkiewicz, 1971) is separately shown 
in Fig. 2 in both the local and physical frames of reference. 
Within a typical element (e)> let the spatial coordinates be 
interpolated as follows: 


= v)z„ =]>]*,( r,„)r, 

/ * l f* 1 

where N, are quadratic “shape” functions associated with the 
element corner, midside and interior nodes. Next, the flow 
variables are interpolated throughout the element in a similar 
fashion. Guided by the Ladyshenskaya-babuska-brezzi com- 
patibility requirements (Carey and Oden, 1986) for the problem 
at hand, the velocity components and pressure are expressed 
as follows: 

= v)V,„ 

/* I ' * I 

K' = Y, .\<f. r,) Y e . „ = £ M* < f, r,)p k 

/ * 1 * “ * 

where M k are the linear shape functions associated with the 
element corner nodes. 

According to the weighted residual method, the error func- 
tions produced by Eqs. (1) through (4) as a result of substituting 
the interpolation expressions, above, are then made orthogonal 
to a special set of weight functions over the finite element 
subdomain. In constructing the latter set of functions, the so- 
called error consistency criterion of Hood and Taylor (1974) 


was implemented, whereby the element shape functions, 
were used in conjunction with the continuity equation. On the 
other hand, quadratic functions which include the element 
shape functions “N" and a set of derived functions “VF” 
were used in conjunction with the momentum equations in 
such a way to ensure full upwinding of the convection terms. 
Of these, the functions ” W ” were previously defined by Hein- 
rich and Zienkiewicz (1977) in terms of the shape functions 
and some upwinding constants which depend on the element 
geometry and local velocity direction. 

With the weight functions now defined, derivation of the 
finite-element equivalent to Eqs. (1) through (4) is straight- 
forward. The process requires linearization of these equations, 
using known values for the velocity components and eddy 
viscosity, and use of Gauss divergence theorem. The final form 
of these equations for the typical element (e) is as follows: 



+ [ f - V rdA p k = <E> £ e rr ( X • V V r )dL 

*1 ^ P dr _ ^(e) 
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h f rdA \p k = b y cf frN,(^ m V V z )dL (13) 
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In these equations, the subscripts and vary from I to 
9, while “fr” varies from 1 to 4. Also the symbol (*) in these 
equations signifies a value that is known from a previous it- 
eration or an initial guess. The global set of equations is achieved 
by assembling Eqs. (11) through (14) among all elements, for 
the current iterative step, and the result is a system of linear 
algebraic equations in the flow nodal variables. 


Results and Discussion 

Two secondary- flow passage configurations, corresponding 
to the same impeller geometry and operating conditions, were 
chosen for comparison. These are shown in Fig. 4, and feature 
a conventional wear-ring and a face seal as part of the sec- 
ondary passage. The pump, which was the focus of rotordy- 
namic testing by Sulzer Bros. (Bolleter et a!., 1989), has the 
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Fig. 4 Sulzer Bros.’ wear ring and lace seal pump configurations (Bol- 
leter et al., 1989) y 
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following characteristics and design-point operating condi- 
tions. 

Impeller tip radius = 17.5 cm 
Impeller speed = 2000 rpm 
Working medium is water at 30°C 
Volumetric flow rate= 130 1/s 
Total head = 68 m 

Reynolds* number (based on the tip speed and ra- 
dius) = 8.02 x 10 6 

In creating the finite-element models in Fig. 4, a total of 
thirteen computational nodes were placed on each cross-flow 
grid line in the seal region. These nodes were closely spaced 
near the walls in anticipation of large velocity gradients there. 
In an early numerical experimentation phase of the study, this 
finite-element grid was proven to provide a good flow reso- 
lution and rule out any significant grid dependency of the 
computed flow field. 

Figure 5 shows a plot of the computed meridional velocity 
component for the conventional wear-ring seal configuration. 
This component, together with the corresponding swirl velocity 
and static pressure (Figs. 6 and 7) constitute the flow’ solution 
corresponding to the current model. Examination of Fig. 5 
reveals that the shroud-to-housing flow is experiencing a pro- 
nounced recirculatory motion in the secondary-passage seg- 
ment leading to the wear-ring seal. This is a result of the 
tendency of the fluid particles adjacent to the shroud to migrate 
radially outwards due to the centrifugal force caused by the 
shroud rotation, on one hand, and the tendency of those par- 





Fig. 7 Contour plot of the nondimensional static pressure for the wear 
ring-seal pump configuration 

tides near the housing to proceed radially inwards as a result 
of the static pressure differential across the passage, on the 
other. 

Contours of the swirl velocity component and static pressure 
associated with the secondary flow field are shown in Figs. 6 
and 7, respectively. The swirl velocity values in Fig. 6 are 
nondimensionalized using the impeller tip speed £/,, and the 
nondimensional pressure (jp) in Fig. 7 is defined as follows: 

__ ip 
P = — 773“ 
pUt 

with p and p l referring to the local and stage-inlet pressures. 
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Fig. 1 0 Contour plot of the nondimensional static pressure for the face- 
seal pump configuration 








Fig. 8 Meridional velocity plot for the face-seal pump configuration 



Fig. 9 Contour plot of the nondimensional swirl velocity lor the face- 
seal pump configuration 


respectively. Also shown in Fig. 7 is a magnified view of the 
wear-ring seal where the major part of the pressure differential, 
across the leakage passage, takes place. 

The meridional and swirl velocity components, along with 
the static pressure obtained for the face-seal pump configu- 
rations are shown in Figs. 8 through 10. Again, the meridional 
velocity vectors, in Fig. 8, indicate a strong recirculatory mo- 
tion in virtually all segments of the secondary passage. Rea- 
soning of this flow behavior was discussed earlier. However, 
a unique, and perhaps peculiar, flow structure is seen to exist 
in the horizontal segment of the leakage passage in Fig. 8, 
where radial shifting of the flow trajectories and vortex break- 
down takes place. It is apparent, however, that the magnitude 
of the pressure gradient in this low-radius region exceeds that 


of the centrifugal force. This seems to weaken the flow recir- 
culation in this leakage-passage segment, confining it to the 
shroud side of the leakage passage. 

There is a rather modest amount of experimental data to 
validate the computed flow field in Figs. 5 through 10. First, 
it was indicated by Sulzer Bros, that the average swirl velocity 
component at the leakage-passage inlet station was measured, 
for the face-seal pump configuration under the above-men- 
tioned operating conditions, to be approximately 0.5 of the 
impeller tip speed (Childs, 1989). With the current finite-ele- 
ment model, the average value of this velocity component was 
computed to be 0.526 of the tip speed. Note that the shroud- 
to-housing swirl velocity profile, w'hich gave rise to this average 
value, is not specified a priori, but is rather part of the finite- 
element flow solution (Fig. 9). This advantage, as indicated 
earlier, is a result of including primary flow segments in the 
computational domain definition (Fig. 4). Further experimen- 
tal observations which are consistent with the computed flow 
field in this study concern the recirculatory pattern of the 
meridional flow (Fig. 8), and were qualitatively reported by 
Guelich et al. (1989). 

Rigorous verification of the seal friction resistance was not 
quite attainable. This was particularly the case for the face- 
seal configuration, which is rather an uncommon seal type. 
Indeed the only “vaguely” similar problem in literature is that 
of a rotating disk in a stationary chamber (e.g.. Daily and 
Nece, 1960) or in open space (e.g., Uzkan and Lipstein, 1986). 
Such studies were focused on the pumping of secondary, often 
cooling, flow over the face of what is typically a bladed-rotor 
disk between the shaft and the disk tip radius. The flow di- 
rection in this case would naturally be radially outwards, in 
contrast to that of the face seal flow (Fig. 4), which is pre- 
dominantly driven by a strong inward pressure decline. Vali- 
dation of the wear-ring seal results, on the other hand, may 
appear possible in light of such experimental correlations as 
those of Yamada (1962). The fact, however, is that Yamada’s 
study did not account for the seal-inlet preswirl velocity which, 
by reference to Fig. 6, is hardly confined to a thin layer adjacent 
to the rotor, nor did the study address the complex flow pattern 
at the end of the shroud-to-housing gap leading to the seal 
segment (Fig. 4). In fact, the combination of these two real 
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flow effects lead to an early, but localized, flow separation 
over the shroud surface at the point marked “S” in Fig. 6, a 
situation that was not encountered in Yamada’s study, where 
only simple isolated annular seals were tested. 

Fortunately, the authors were recently able to validate the 
current computational model, including the turbulence closure, 
for a straight annular seal using Yamada’s experimental data 
(Baskharone and Hensel, 1991). The seal under consideration 
had a clearance/length ratio of 0.034, and a Reynolds’ number 
(based on the clearance width and inlet through-flow velocity) 
of 13,280. The computed seal resistance for this seal was re- 
ported to be in excellent agreement with Yamada’s experi- 
mental data. Furthermore, the agreement between the 
computed profiles of through-flow and tangential velocity 
components and those of Morrison et al. (1988) was equally 
favorable. 

Performance assessment of the two secondary- flew passages 
(Fig. 4), in light of the numerical results, involves their effec- 
tiveness as leakage suppressants and swirl dissipators. In order 
to determine the leakage control capacity of each passage, the 
mass flux was integrated at the passage inlet station. The results 
indicated that the leakage-flow rates in the w r ear-ring and face- 
seal pump configurations were 0.0011 and 0.0038 m 3 /s, re- 
spectively. These represent 0.85 and 2.92 percent of the primary 
flow rate, and illustrate the relative superiority of the wear- 
ring seal configuration, as a leakage-control device over the 
face-seal alternative. This, in part, is due to the highly favorable 
streamwise pressure gradient across the face seal as a result of 
the substantial radius diminishmem, between the seal inlet and 
exit station, w r hich lessens the boundary layer build-up over 
the solid walls and eliminates the likelihood of any flow sep- 
aration. On the other hand, the tendency of the fluid particles 
to migrate radially outwards near the constant-radius shroud 
surface, in the case of the wear-ring seal, added to the complex 
flow structure in the passage leading to the seal, create an 
environment for an early flow separation over the shroud sur- 
face (location “S” in Fig. 6). Despite the rapid flow- reattach- 
ment, in this case, the recirculation region, following the 
separation point, has the effect of enhancing the seal frictional 
resistance and, therefore, the sealing effectiveness of the wear- 
ring seal configuration by comparison. Another contributing 
factor behind the wear-ring effectiveness, as a leakage-control 
device, is the existence of a substantial recirculation zone at 
the seal-exit station (Fig. 5), with practically no comparable 
seal-exit flow behavior in the case of the face-seal configuration 
(Fig. 8). As for the swirl velocity dissipation across the sec- 
ondary passage, examination of Figs. 6 and 9 reveals that the 
“kink” in the secondary passage of the face-seal pump con- 
figuration (Fig. 9), is causing a sudden and measurable re- 
duction in the swirl velocity in the axial passage segment, with 
no equivalent swirl velocity reduction mechanism in the wear- 
ring seal pump configuration. More importantly, the average 
swirl velocity component at the secondary passage inlet station 
(designated x-x in Fig. 4) was found to be as low as 0.526 of 
the impeller tip speed for the face-seal pump configuration, 
as opposed to 0.812 for its counterpart. Worth noting is the 
fact that the choice of station x-x in Fig. 4 is consistent with 
that of Childs (1989), as a meaningful parameter in interpreting 
the rotordynamic behavior of the same pump configurations 
in Fig. 4. Although these swirl characteristics of (he face-seal 
pump configuration are not quantitatively convertible into a 
rotordynamic stability-related factor, it is established (Childs 
et al., 1991) that swirl suppression at the secondary passage 
inlet station is among the most effective tools for shrouded- 
impeller rotordynamic stability enhancement. This, in view of 
the current results, would imply that the face-seal pump con- 
figuration provides a less destabilizing effect. The conclusion 
here is consistent with the experimental findings of Bolleter et 
al. (1989), who reported comparable direct-damping coeffi- 
cients for the fluid/shroud interaction system, but a signifi- 


cantly lower cross-coupled stiffness coefficient for the face- 
seal pump configuration. Since unstable operation of the im- 
peller (in the form of a whirling motion) is triggered by low 
direct damping and/or high cross-coupled stiffness, the ex- 
perimental measurements suggest that the face-seal pump con- 
figuration results in a more stable impeller operation, which 
is what the numerical results seem to imply. 


Concluding Remarks 

The finite-element model of the pump secondary flow in 
this paper is a versatile tool for predicting the leakage-suppres- 
sion and rotordynamic-stability characteristics of viable sec- 
ondary-passage designs. The model accuracy stems from the 
manner in which the computational domain is defined to permit 
the typically strong primary/secondary flow interaction. Sim- 
ulation of the flow turbulence and recognition of the inertia 
domination in formulating the problem makes the model ap- 
plicable to a wide range of real-life operating conditions. The 
model was applied to two representative configurations of seal- 
containing secondary passages in an existing centrifugal pump 
stage for code verification and seal-performance assessment 
purposes. The results are consistent with experimental meas- 
urements concerning the same pump, and documented obser- 
vations in similar secondary-passage configurations. The 
computed flow field in the current study constitutes the so- 
called “zeroth-order” flow solution in existing perturbation 
models of the fluid-induced vibration of shrouded pump im- 
pellers. This explains the desire for an accurate flow field in 
the secondary passage, since the output of such perturbation 
models would naturally be a strong function of the flow’ be- 
havior prior to the impeller excitation. 
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Introduction 


Frequent mechanical failure of turbomachinery rotating components, due to fluid induced 
vibration, has been reported over the past two decades. The most common instability mechanism 
in this case is one where the impeller whirls around the housing centerline at a subsynchronous whirl 
frequency (e.g. Childs and Moyer, 1985). Serious efforts are now being devoted to the development 
of reliable predective tools for capturing the dynamics of fluid-encompassed rotors and testing ways 
to enhance the fluid restoring forces. 

Rotordynamic forces on shrouded pump impellers have particularly been the subject of an 
increasing number of experimental and numerical investigations. Of the experimental studies, 
those by Bolleter et al. (1989) and Guinzburg (1992) are perhaps the most comprehensive. In both 
cases, the fluid-exerted forces were reported to have a more or less parabolic dependency on the 
impeller whirl frequency. This dependency was approximated with a least-square fit to obtain the 
direct and cross-coupled stiffness, damping and inertia coefficients of the fluid/shroud system. 

Perhaps the simplest, and most popular, numerical tool in this area is the bulk -flow model 
devised by Childs (1983). This model was later upgraded (Childs, 1989) and used to compute 
the rotordynamic coefficients of two leakage-passage configurations, both corresponding to the 
same impeller geometry and operating conditions. Results of this study were compared to the 
experimental data by Bolleter et al. (1989). More detailed models of the fluid/rotor interaction 
problem have also been reported by Tam et al. (1988) as well as Dietzen and Nordmann (1988). 
However, attention in these two studies was focused on simple annular seals and journal bearings. 

The current study is the first such attempt, known to the authors, where a full-scale unsplit 
secondary passage of a pump stage is under investigation within a framework which tolerates such 
real flow effects as separation and recirculation. The study is an extension of the model previously 


published by Baskharone and Hensel (1991a) for a simple untapered annular seal. The idea of uti- 
lizing segments of the primary flow passage in defining the computational domain (Fig. 1) was first 
outlined (also by Baskharone and Hensel, 1989) as a means of facilitating the primary /secondary 
flow interaction effects. However, the published analysis, then, merely presented the idea, but 
treated the secondary flow as strictly laminar and, as a result, had no appreciable engineering 
value. The current study is, therefore, one which expands previously devised ideas, by the authors, 
and applies the outcome to a traditionally complex problem for which an equally detailed solution 
is simply non-existing. 

Analysis 

Fig. 1 shows a schematic of a hydraulic pump stage, with a face seal as a leakage-control device. 
Shown on the same Figure is the computational domain for solving the centered-rotor (zeroth- 
order) flow-governing equations. In addition to the secondary shroud-to-housing flow passage, 
the selected domain also includes two primary-flow segments, which are naturally connected to 
the secondary passage near the impeller inlet and discharge stations. The objective here is to 
alleviate the need for specifying, assuming, or iteratively obtaining the secondary-passage inlet 
flow conditions. Of these common choices, Childs, for instance, has reportedly utilized the iterative 
approach, whereby the secondary passage was isolated, and the flow solution iteratively obtained 
under the condition that the streamwise pressure differential, across the passage, is identical to 
that across the impeller (Childs, 1989). Such approach is understandably appropriate, considering 
that the secondary flow in this case is treated as one layer with no cross-flow gradients taken into 
account. However, this and any other equivalent simplifications would fall between insufficient and 
inadequate in the current detailed flow model where inlet and exit velocity profiles will also have to 
be prescribed, should the passage be isolated. Contouring the computational domain in the manner 
seen in Fig. 1 has the effect of converting the the secondary passage inlet and exit flow conditions 
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into internal variables to be part of the numerical solution. Note that the passage inlet and exit 
swirl profiles, which have a predominant effect on the impeller stability, are included among these 
variables. 

Fig. 2 shows the finite-element discretization model which was created for the centered-rotor 
operation mode. The model is composed of non— overlapping nine— noded bi— quadratic elements 
of the Lagrangian type (Zienkiewich, 1971). These elements are generally curve-sided, a feature 
which had a highly favorable effect in matching the domain boundaries. Intensity of the finite 
elements are, by reference to Fig. 2, much higher in regions where pronounced gradients of the flow 
thermophysical properties are anticipated. 

Flow— governing Equations. In order to eliminate the flow field time-dependency, that is 
imparted by the impeller whirling motion (Fig. 3), the flow conservation laws are cast in a frame 
of reference that is attached to the rotor, being the impeller-shroud assembly, and whirls with it. 
The time-avaraged flow equations, now expressed in terms of relative velocity components, include 
such terms as the Coriolis and centripetal acceleration components. These equations, as well as 
their finite-element equivalent, were previously published by Baskharone and Hensel (1991b). The 
published model also contained a modified version of Baldwin and Lomax turbulence closure (1978). 

Boundary Conditions. Referring to the flow-permeable boundary segments in Fig. 1, the 
boundary conditions over these segments are as follows: 

i) Stage Inlet Station: This is the boundary segment a-b in Fig. 1, which is located sufficiently 
far upstream from the impeller. Fully developed flow is assumed at this location, giving rise to the 
following boundary condition: 


dV r _ dV e _ dV z 

dz dz dz 
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In addition, the stage-inlet static pressure is specified at the node midway between the endwalls 
on this station. 

U) Impeller Inlet and Exit Stations: These are labled c-d and e-f in Fig. 1. Fixed profiles of 
the velocity components, corresponding to the stage operating conditions, are imposed over these 
boundary segments. Note that the operating conditions here involve the primary impeller passage, 
and do not include the secondary mass flow rate. 

m) Stage Exit Station: The flow behavior at this station (designated g-h in Fig. 1) is viewed 
as predominantly confined to satisfying the mass and angular momentum conservation equations 
in a global sense. In their derivative forms, These can be expressed as follows: 


OVr 

dr 



and 


dV e 

dr 



These two boundary conditions are linear and are, therefore, introduced non-iteratively in the 
numerical solution process. Moreover, a zero normal derivative of V z is imposed over this station, 
and the stage-exit static pressure is specified at the computational node midway between the 
endwalls on this station. 


As for the solid boundary segments in Fig. 1, namely those of the housing and shroud as well 
as the hub surface segment b-d, the no-slip boundary condition applies as follows: 


V r = 0, V z = 0 and V e = C 

where C is equal to (fir) and zero for rotating and non-rotating boundary segments, respectively. 

Perturbation Model. As the impeller enters the whirling motion depicted in Fig. 3, the 
flow properties including, in particular, the shroud pressure distribution, assume non— axisymmetric 
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patterns. The current perturbation model perceives such deviation, from that of the centered- 
impeller, as the result of infinitesimally small distortions in the finite-element grid due to the 
eccentricity V J of the impeller axis of rotation (Fig. 3). Expansion of the flow-governing equations, 
in their finite-element form, as functions of V 5 is then carried out and, in the end, gives rise to 
the differential changes of the flow variables. Of these differential changes, those associated with 
the shroud pressure distribution are then isolated and integrated over the entire shroud surface. 
The result, in this case, is the rate, with respect to e at which the fluid reaction forces (restoring 
or aggravating) are exerted on the shroud. The fluid/shroud interaction forces are then analyzed 
at various whirl frequencies to determine the stiffness, damping and inertia coefficients with which 
the fluid contributes to the shroud whirling motion. 

Aside from the different geometry and boundary conditions in the current study, the perturba- 
tion model just described is virtually identical to that by Baskharone and Hensel (1991) where the 
much simpler case of a straight annular seal was analyzed. Note that both of virtual eccentricity 
“e” and the whirl frequency “fl”, do affect the differential changes in boundary conditions in the 
process of expanding the finite-element equations. 


Results and Discussion 

A typical pump stage, designed and tested by Sulzer Bros. Inc. (Switzerland), was chosen 
for the current study. This is schematically shown in Fig. 1, and has been previously selected 
for rotordynamic analysis by Childs (1989). The stage dimensions and operating conditions were 
reported by Bolleter et al. (1989) as follows: 

Impeller tip radius = 17.5 cm 


Impeller speed = 2000 rpm 



Working medium is water at 30 °C 

Volumetric flow rate = 130 Ijs 

Total head = 68 m Reynolds 1 number (based on the tip speed and radius) — 8.02 x 10 6 

Figures 4, 5 and 6 present the most significant features of the centered-impeller flow structure 
in the secondary passage. The complexity of this flow field is most apparent in Fig. 4, which is 
a vector plot of the meridional velocity component in the shroud-to-housing passage. The flow 
pattern in this Figure is characterized by massive separation and recirculation zones in virtually 
all segments of the passage. This characteristic is the result of the flow tendency to migrate 
radially outwards near the “spinning” shroud (due to the locally high centrifugal force and the flow 
viscosity), which is opposed by the tendency to proceed radially inwards near the housing (due to 
the static pressure differential across the flow passage). Such flow behavior is generally consistent 
with experimental observation (e.g. Guelich et al., 1987), but was never part of a rotordynamic 
analysis to the authors’ best knowledge. Equally important are the swirl velocity profiles in Figs. 
5 and 6 at the leakage-passage and the seal flow inlet stations (Figs. 5 and 6, respectively). 
This is due to the well-established fact that the flow swirl in the leakage passage acts as a major 
destabilizing factor. Examination of the cross-flow swirl velocity distributions in Figs. 5 and 6 
reveals significantly large boundary layer thicknesses at the solid walls, with a rather narrow “core 
flow” region as the flow passage contracts just upstream of the face seal (Fig. 6). The relatively 
small value of the swirl velocity in Fig. 6 is a natural outcome of the “no-slip” condition which 
prevails at the shroud surface which is relatively small at this radial location by comparison. 

Prior to the full-scale execution of our perturbation model, a grid-dependency study was 
conducted. The independent variable here was the number of computational planes in the circum- 
ferential direction (Fig. 3), whereby the fluid-exerted forces at an arbitrarily fixed whirl-frequency 
ratio (fl/u;) of 0.3, were monitored. Results of this numerical experimentation phase are shown in 


Fig. 7 in the form of tangential and radial components of the net shroud force. The Figure clearly 
shows that the force trends begin to level out for a circumferential plane count of 9. As a result, 
the impeller rotordynamic analysis was conducted in its entirety using this value in a compromise 
between the precision of the numerical results, on one hand, and the consumption of computational 
resources, on the other. 

Results of the perturbation analysis, as fluid-exerted forces on the shroud, are shown in Fig. 
8 for a range of whirl frequency ratio (H/u;) between -1.25 and +1.25. The fluid-exerted forces in 
this Figure, as well as Fig. 8, are non-dimensionalized in the same manner proposed by Bolleter 
et al. (1989), namely: 


0F r 

d€ 

7?r t 2 pb t uf 2 


QF t 

p 9e 

£ Ttr t 2 pb t uj 2 

where € is the orbit radius of the impeller axis (Fig. 3), r t is the impeller tip radius, b t is the 
impeller-tip endwall spacing, p is the fluid density, and u; is the running speed. Reproduced on the 
same Figure are the experimental measurements of Bolleter et al. (1989), and the forces computed 
by Childs (1989) using a bulk-flow perturbation model. Of these, Childs’ results are shown for 
pre-imposed inlet swirl values of 0.5 and 0.6 of the impeller tip speed. The reader is reminded that 
this ratio, in the current computational procedure, is treated as unknown to be produced in the 
solution of the zeroth-order flow field, and that our computed average value is, by reference to Fig. 
6, 0.526 of the impeller tip speed. As seen in Fig. 8, the bulk-flow model produces rather sharp 
force fluctuations, initially thought of a resonance-like, for the elevated inlet swirl ratio of 0.6, 
within a narrow range of whirl-frequency ratio of 0.3. These fluctuations are neither confirmed by 


the current results or the experimental data. The fact that such drastic change in the force trends 
corresponds to modest changes in the inlet swirl ratio, by comparison to the smooth flow trends 
obtained for a swirl ratio of 0..5 (Fig. 8), was equally surprising. However, the point should be made 
that the seal-inlet swirl velocity in Childs’ study was significantly higher than (roughly double) the 
computed value in the current study (Fig. 6). Knowing that the seal-inlet swirl velocity deviation, 
between the two numerical studies, grows even wider as the leakage-passage inlet swirl velocity is 
elevated, it would perhaps be natural to expect the type of disagreement between the two sets of 
results in Fig. 8. 

Next, the direct and cross-coupled stiffness, damping and inertia (added mass) coefficients 
of the fluid/shroud interaction were computed. The computational procedure here is essentially 
equivalent to that previously outlined by Baskharone and Hensel (1991a) for an annular seal prob- 
lem. The calculated rotordynamic coefficients are shown in Table 1 versus the experimental results 
of Bolleter et al. (1989) and the numerical results of the bulk-flow model (Childs, 1989). Of these 
the results by Childs correspond to an externally-prescribed leakage-passage inlet swirl ratio of 0.5 
which, by reference to Fig. 8, gives rise to smooth trends of the fluid reaction force components. 
Referring to Table 1, it seems that the bulk -flow model significantly underestimates the direct 
stiffness and cross-coupled inertia coefficients, particularly the latter. Otherwise, the rotordynamic 
coefficients produced by the two perturbation analyses appear to be of comparable magnitudes. 

In appraising the current numerical results in Table 1, the fact should be emphasized that a 
significantly better agreement with the experimental data was not initially anticipated. Despite our 
success in facilitating (to a reasonable extent) the primary/ secondary flow interaction within the 
pump stage (Fig. 1), it is the intolerance of our perturbation model to the force contributions of the 
impeller primary flow which led to such modest expectations in the first place. Nevertheless, the 
fact remains that our finite-element -based perturbation model, including the relatively rigorous 
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flow analysis on which it is founded, is by no means less accurate than any existing perturbation 
model known to the authors. 

Concluding Remarks 

A perturbation model, originally devised for simple annular seals, has been upgraded and 
applied to the problem of fluid/shroud interaction in a typical pump stage. The computed rotordy- 
namic coefficients are compared to existing experimental and numerical data. The current results 
seem closer to the experimental values for those coefficients where a substantial disagreement, 
with the existing numerical data, occurs. A primary source of discrapencies between the newly- 
computed sets of rotordynamic parameters and their numerical counterparts (Childs, 1989) was 
identified to be the seal-inlet swirl velocity magnitude. In the current analysis, this key variable is 
part of the zeroth-order (centered-impeller) flow solution and, as such, is a function of not only the 
pump operating conditions but also the manner in which the leakage passage is naturally connected 
to the seal segment of the secondary passage. The authors believe that the value of the current 
perturbation model would be enhanced by integrating the impeller primary flow passage into the 
computational domain. Such upgrade, which is major in nature, is currently under consideration. 
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Fig. 5 Distribution of the non-dimensional swirl velocity component 
at the lealage-passage inlet station 



Fig. 6 Distribution of the non-dimensional swirl velocity component at the face-seal inlet station 
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Fig. 7 Dependency of the fluid-exerted forces on the grid resolution in the circumferential direction 
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Fig. 8 Comparison of the shroud forces with experimental and numerical data 
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ABSTRACT 

A detailed fmite-element-based perturbation analysis, initially developed for straight annular 
seals, is extended to the labyrinth se?! rotordynamics problem. This particular seal category is 
frequently used in turbopumps, such as those of the Space Shuttle Main Engine, for its effectiveness 
as a leakage control device. The same seal configuration is also known to produce a net destabilizing 
force, which stems from the fluid/rotor interaction. The general objective here is to assess the 
effect of a design parameter, namely the number of tooth-to— tooth chambers, on the seal relative 
stability. To this end, two representative cases of teeth-on-rotor labyrinth seals are selected for the 
purpose of comparison. The chosen seals share the same tooth-to-tooth chamber geometry, but are 
composed of one and five such chambers. The centered-rotor flow solution corresponding to each 
seal version is first obtained, and relevant flow characteristics in each case outlined. This “zeroth— 
order” flow solution is next used as input to the perturbation model, with the outcome being the 
fluid-exerted forces on the rotor as the latter whirls around the housing centerline. These forces are 
then analyzed to determine the direct and cross-coupled stiffness, damping and inertia coefficients 
of the fluid/rotor interaction. The computed rotordynamic coefficients are compared to an existing 
set of experimental measurements concerning the same two seal configurations. The final outcome 
of the study is consistent with the general belief that labyrinth seals impart a rotordynamically 
destabilizing effect. In terms of relative stability, however, the five-chamber seal configuration offers 
less of this destabilizing effect, by comparison. 
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NOMENCLATURE 




C ~ direct damping coefficient of the fluid/rotor interaction 
c = cross-coupled damping coefficient 
F r = radial component of the fluid-induced force 
F r ~ = dF r /de 

F$ = tangential component of the fluid-induced force 
F 0 m = dF 0 /d€ 
f — whirl frequency ratio 
h = seal clearance 

K = Direct stiffness coefficient of the fluid/rotor interaction 
k — cross-coupled stiffness coefficient 
p = static pressure 
Pi = seal-inlet static pressure 
r = radius 

r* — inner radius of the labyrinth seal 
V z — axial velocity component 
V$ = tangential velocity component 
e = lateral eccentricity of the rotor axis 
p — density 
f2 = whirl frequency 
w — rotor operating speed 
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INTRODUCTION 


Labyrinth seals are generally regarded as efficient leakage-control devices in turbomachinery 
applications (Fig. 1). Since leakage flow is considered on of the primary loss mechanisms, and 
subsequently performance degradation, in turbomachines, much of the numerical studies on seals 
have been focused on this particular seal category. Among those are the computational leakage 
models by Stoff (1980), Rhode et al. (1984) for incompressible flow applications. Compressible-flow 
labyrinth seal models, for gas turbine applications, have also been devised le.g. Wittig et al. (1987) 
and, more recently, Rhode and Hibbs (1992)]. 

Experimental measurements and flow visualization studies in labyrinth seals have also been re- 
ported in the literature. Examples of these include the LDA measurements by Wittig et al. (1987), 
the flow visualization study by Iwatsubo and Kawai (1984), and the discharge coefficient measure- 
ments by Wittig et al. (1987). The complexity of the flow structure in labyrinth seals has also 
encouraged researchers to devise, and continually update, empirical correlations (e.g. Zimmermann 
and Wolff, 1987). 

Despite the performance improvements attained through the utilization of labyrinth seals, se- 
rious fluid-induced vibrations have been attributed to them by many rotordynami cists. Among 
those, Alford (1965) was the first to point out the destabilizing effects of labyrinth seals. Iwat- 
subo and Kawai (1984) later reported significantly poor direct damping coefficients for a hydraulic 
labyrinth seal under a matrix of different operating conditions. In a comparative study involving 
several seal categories, Childs and Elrod (1988) also concluded that the flow swirl at the labyrinth 
seal inlet station would make a bad rotordynamic situation even worse, and that an effective swirl 
brake would be required in this case. 

Recent theoretical developments in this area have generally been either too simplified to ac- 
count for real flow effects in labyrinth seals, or so tailored as to handle simple geometries and/or 
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uniform lateral eccentricity of the whirling rotor. Under the first category are the “one-volume” 
models by Iwatsubo (1980) and that by Childs and Scharrer (1984). To the author’s knowledge, 
the first attempt to apply rigorous computational fluid dynamics tools, within the framework of 
a perturbation approach, was that by Xordmann and Weiser (1988). This finite difference-based 
model is based on the transformation of the entire rotor-to-housing computational domain to a 
fictitious frame of reference whereby the rotor eccentricity (assumed uniform along the seal axis) 
is eliminated. Despite the apparent generality and complex nature of this approach, the method 
is practically limited to rectangular tooth-to-tooth chambers, and is concepuallv incapable of ad- 
dressing rotor excitations in the form of conical whirl, which would be of serious consequences for 
long seals with appreciable streamwise pressure differential. 

The current perturbation model is a versatile predictive tool, which is based on what we 
generically termed the “virtually” deformable finite-element concept. Theoretical details of this 
model were documented by Baskharone and Hensel (1991a), and a sample case of a straight annular 
seal successfully analyzed (Baskharone and Hensel, 1991b). Versatility of this model was also 
demonstrated by analyzing, as interrelated, the cylindrical and conical whirl of a straight annular 
seal rotor (Baskharone and Hensel, 1991c). The present study illustrates another versatility aspect 
of the perturbation model, namely that of the seal geometry arbitrariness. The study addresses 
both the traditional difficulty of analyzing the flow past a sequence of high aspect-ratio cavities, 
normally found in labyrinth seals, as well as the task of classifying, as restoring or aggravating, the 
fluid reaction forces due to excitations of the rotor axis (Fig. 2). 


COMPUTATIONAL DEVELOPMENT 

12 Zeroth— Order Flow Field. The perturbation procedure is initiated by computing this 

“zeroth-order” flow field, which is clearly axisymmetric. Details of the flow-governing equations, 
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turbulence closure, and finite-element formulation of this relatively simple problem were all dis- 
cussed by Baskharone and Hensel (1991b), where a test case involving a straight annular seal was 
comprehensively covered. Comparison of the computed through— flow and swirl \elocit\ profiles in 
the rotor-to-housing passage with the experimental data previously reported by Morrison et al. 
(1991), was also presented, and the agreement between the two sets of data was encouraging. 

Perturbation Analysis. The centered-rotor flow solution, referenced above, was used 
in the current study as input to the perturbation analysis phase of the computational procedure. 
Devised by Baskharone and Hensel (1991a), the perturbation model here is based on what was 
generically termed the “virtually” deformable finite element concept, where the perturbed flow 
equations emerge from expansion of the finite— element equations in terms of the rotor eccentricity 
e. Since no circumferential pattern is pre-imposed on the perturbations of the flow properties, the 
problem at this point shifts to the fully three-dimensional type, as the rotor eccentricity destroys 
the flow axisymmetry (Fig. 2). Simultaneously, the flow problem is cast in a rotating-translating 
frame of reference, which is attached to the whirling rotor at all times (Fig. 2) in order to eliminate 
the time dependency of the flow field in this case. 


RESULTS AND DISCUSSION 

Two labyrinth seals, composed of the same tooth-to-tooth chamber geometry, but different in 
the number of chambers were selected for this study. These one and five-chamber hydraulic seals 
were the subject of a flow visualization study and rotordynamic testing by Iwatsubo and Kawai 
(1984). In both cases, a seal-inlet pressure and rotor speed of 147 KPa and 3.7 Hz were arbitrarily 
chosen from the experimental matrix of operating conditions. The finite-element model created 


for the two seal configurations are shown in Figs. 3 and 4, respectively. The discretization unit 
in these two Figures is a bi-quadratic curve-sided element of the Lagrangian type. The process 
of generating finite-element grids, such as these, was largely automated and the level of near-wall 
refinement in the preprocessing segment of the computational procedure was made totally arbitrary. 
This capability made it possible to conduct a preliminary numerical “experimentation” study of 
grid optimization for the centered-rotor (axisymmetric) flow solution. The finite-element grids in 
Fig. 3 and 4 are the outcome of this early study, where a careful compromise was made between the 
numerical accuracy, on one hand, and the consumption of computational resources, on the other. 

Centered-Rotor Flow Field. Figures 5 through 7 show the centered-rotor flow solution 
for the one-chamber seal configuration. This solution constitutes the zeroth-order flow field in the 
current perturbation analysis. First, the meridional flow behavior is shown in Fig. 5, with enlarged 
segments of the computational domain to clarify the flow recirculation wdthin the tooth-to-tooth 
chamber and the vorticity breakdown in the dump region downstream from the chamber. The flow 
pattern in this Figure was qualitatively compared to that resulting from a flow visualization study 
by Iwatsubo and Kawai (1984), and the major features of the computed flow pattern were found 
to be consistent with the experimental findings. 

The swirl velocity contours, throughout the computational domain, are shown in Fig. 6. 
The swirl velocity (V$) in this Figure is non-dimensionalized using the average inlet through-flow 
velocity component (V,;). Note that The Figure shows no inlet preswirl, which is consistent with 
Iwatsubo’s seal-inlet conditions. Also note the swirl-velocity boundary layer profile development 
over the rotor surface. As seen in the Figure, the swirl velocity near the housing surface appears in 
the tooth-to-tooth chamber region, and steadily declines, over the housing, away from the chamber, 
up to the location where the flow recirculation zone, downstream from the second tooth (Fig. 5) 


Fig. 7 shows a contour plot of the static pressure (p), with an enlarged view of the tooth-to- 
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tooth seal chamber. The non-dimensional pressure (p) in this Figure is defined as follows: 


P = 


(p - pd 

pV :i 2 


where p and p x are the local and seal-inlet static pressures, respectively, p is the fluid density, and 
V zi is the average through-flow velocity at the seal inlet station. Note that the upstream side of 
the second tooth is exposed to a rather high pressure near the housing, despite the fact that the 
center of the low-pressure island is closer to the that tooth. 


The centered-rotor flow field corresponding to the five-chamber seal configuration is shown in 
Figs. 8 through 10, with particular emphasis on the flow structure in the first and last chambers. 
Examination of Fig. 8 reveals that the center of the cavity recirculatory motion gradually moves 
from the cavity center towards the tooth pressure side as the flow progresses from one cavity to the 
next. The most significant difference between the swirl velocity distribution in Fig. 9 and that of the 
single-chamber seal (Fig. 6) is the continuous rise in the magnitude of this velocity component in 
the through-flow direction, by comparison. This characteristic simply means that the five-chamber 
seal gives rise to “elongated” flow trajectories over the last few chambers. One would expect, as 
a result, a substantial total pressure differential across these chambers, which is consistent with 
the known fact that the leakage-suppression capability of labyrinth seals is drastically enhanced by 
increasing the number of tooth-to-tooth chambers. As for the static pressure distribution in Fig. 
10, it seems that the largest pressure increments occur across the first and last chambers. 


Grid Dependency, 

As would naturally be anticipated, the perturbed flow resolution is dependent on the grid 
C resolution which, as the domain three-dimensionality prevails, translates into the number of com- 

putational stations in the circumferential direction (Fig. 2). Investigation of this dependency was 
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carried out by varying this number and repeatedly computing the fluid-exerted forces on the rotor 
for an arbitrarily chosen whirl frequency "fi” that is identical to the operating speed “u;” . Results 
of this preliminary step are shown in Fig. 11, where the fluid force components F/ and are ob- 
tained by integrating the rotor-surface pressure forces, upon resolution in the radial and tangential 
directions, over all of the contributing finite element boundaries. 

The force trends in Fig. 11 suggest that the radiai force are more sensitive to the tangential 
grid resolution by comparison. The Figure also shows that changes in the force magnitudes becomes 
acceptably small as the tangential computational-station count approaches 11. Given the fact that 
F r itself varies by as low as 2.7 % by increasing the number of these computational stations from 
10 to 11, and in an attempt to maintain a “manageable” CPU time consumption, the number was 
fixed at 11 thereafter. 

Fluid-Induced Forces and Rotordynamic Coefficients, The trends of the fluid- 
exerted forces are shown in Fig. 12 for the five-chamber seal configuration over a range of the 
rotor whirl frequency between -200% and 200%, with the negative values representing backward 
whirl. The reassuring feature of the radial and tangential force trends in this Figure is that each 
of them is more or less parabolic, a commonly known characteristics of seals in general. Of the two 
components in Fig. 12, the tangential force is exlusively responsible for aggravating or suppressing 
the rotor whirling motion. With this in mind, the impression one would have, by examining Fig. 12, 
is that of rotordynamic instability in the range of positive (forward) rotor whirl frequency between 
zero and 20% of the rotor speed, judging by the positive tangential force in this range. 

The computed fluid-exerted forces were then used to compute the rotordynamic coefficients of 
the fluid/rotor interaction, in the manner previously outlined by Baskharone and Hensel (1991a). 



The results of this computational step are contained in Tables 1 and 2 for the one- and five- 
chamber seal configurations, respectively. Reproduced in these tables are also the experimental 
measurements by Iwatsubo and Kawai (1984) for comparison. Perhaps the most alarming values 
in these two tables is the poor value of the direct damping coefficient “C” in both cases. In fact, 
the negative value of this coefficient, for the single-chamber seal configuration, implies a highly 
destabilizing effect of this seal, particularly in view of the positive (also destabilizing) value of the 
cross-coupled stiffness coefficient “k” in this case. 

A useful parameter in quantifying the relative rotordynamic stability of general fluid-encompassed 
rotors is the so-called whirl frequency ratio “f”, where: 



This is the ratio of the destabilizing-to-stabilizing forces acting on the rotor. A given seal would 
have a stabilizing influence for a whirl frequency ratio “f” below 1.0 (Scharrer, 1988). Application of 
this criterion to the five-chamber seal configuration leads to the conclusion that this particular seal 
imparts rotordynamic stability to the shaft (in a global sense), for it gives rise to a whirl frequency 
ratio of approximately 0.22. However, the “relative stability” criterion, above, is inapplicable to the 
single— chamber seal configuration. This is due to the negative value of the direct damping coefficient 
in this case. Such characteristic simply renders the seal unstable regardless of the magnitude of the 
cross-coupled stiffness coefficient, should the latter be positive. 


CONCLUDING REMARKS 


The objective of this study was twofold; to demonstrate the applicability of a categorically 
new perturbation model to the traditionally complex problem of labyrinth— seal rotordynamics, and 



to assess the stability effects of multiple-chamber seal configurations. Of these, the latter task 
is largely design-related, and is currently receiving an increasing amount of attention, owing to 
documented fluid-induced vibration problems in the Space Shuttle Main Engine turbopumps, for 
which labyrinth seals are primarily blamed. An intermediate outcome of the current stud\, namelv 
the centered-rotor flow field, is equally suggestive from a design standpoint. This is true in the 
sense that knowledge of the flow pattern, for this operation mode, provides a means of predicting 
the seal effectiveness as a leakage-control device. 
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Fig. 2 Distortion of the rotor-to-housing finite element assembly as a result of the rotor eccentricity 
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Fig. 3 Finite-element model of the single-chamber seal configuration 
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Fig. 4 Finite-element model of the Ave-chamber seal conBguration 
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Fig. 7 Contour plot of the static pressure for the centered-rotor operation mode 
of the single-chamber seal configuration 




Fig. 8 Vector plot of the meridional velocity component for the certered-rotor operation mode 
of the five-chamber seal configuration 
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Fig. 11 Dependency of the fluid-exerted forces on the grid resolution in the circumferential direction 
for the case of synchronous whirl (Cl/ui = 10) 
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Fig. 12 Fluid-exerted forces on the whirling rotor 
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Table 1 Comparison of the rotordynamic coefficients with experimental data 
for the single- chamber seal configuration 
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Table 2 Comparison of the rotordynamic coefficients with experimental data 
for the five- chamber seal configuration 












Perturbed Flow Structure in an Annular Seal 
Due to Synchronous Whirl “ 

E,A. Baskharone ** 

Department of Mechanical Engineering 
Texas A&M University 
College Station, Texas 77843 

This paper provides a thorough examination of the flow field resulting from synchronous whirl 
of an eccentric rotor in an annular seal under typical operating conditions, A new finite-element- 
based perturbation model is employed in the analysis, whereby perturbations in the flow ther- 
mophysical properties are attributed to virtual distortions in the rotor-to-housing finite element 
assembly. The numerical results are compared to a recent set of experimental data for a hydraulic 
seal with typical geometrical configurations and a synchronously whirling rotor, Despite the com- 
mon perception that perturbation analyses are categorically confined to small rotor eccentricities, 
good agreement between the computed flow field and the experimental data is obtained for an 
eccentricity/ clearance ratio of 50%. The agreement between the two sets of data is notably better 
at axial locations where the real-rig flow admission losses have diminished, and up to the seal dis- 
charge station. This attests to the accuracy of this untraditional and highly versatile perturbation 
model in predicting the rotordynamic characteristics of this and a wide variety of conceptually 
similar fluid/rotor interaction problems. 

Nomenclature 

c = nominal seed clearance 
d = rotor diameter 
e = eccentricity ratio (e/c) 

* Thi» research was funded by NASA-Marshall Space Flight Center (Huntsville, Alabama), Contract No. NAS8-37821, 
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F r w = radial component of the rotor force perturbation 
Fq m ^tangential component of the rotor force perturbation 
L = seal length 
p — static pressure 
p in = seal inlet pressure 
Re = Reynolds number = 2 cVi n ju 
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T a — Taylor number = 

Ui = rotor surface velocity = uri 
V in — seal inlet through-flow velocity 
V r — radial velocity component 
V z = axial velocity component 
V$ = tangential velocity component 
z = distance along the seed axis 
e = rotor eccentricity 
v — kinematic viscosity coefficient 
p — fluid density 
fl = whirl frequency 
u = rotor spinning speed 
()* implies the operator ^ 
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Introduction 


An unstable operation mode of annular seals is known as cylindrical whirl (Fig. 1), in which 
the rotor axis precesses around the housing centerline at a finite frequency (fi). Of the possible 
vibration mechanisms in this case is synchronous whirl, where the whirl frequency is identical to 
the rotor speed, and is primarily caused by mass imbalance. In high speed pump applications, 
the hydrodynamic forces associated with the rotor whirl can be a major contributor to the system 
stability, and may indeed lead to destructive consequences. 

Existing computational tools for predicting the fluid/rotor interaction forces vary in complexity 
from the simple and widely used bulk-flow model (Childs, 1983) to more detailed finite-difference 
analyses (e.g. Dietzen et al., 1987, and Tam et al., 1988). These have progressively embraced 
more of the flow details in the rotor-to-housing passage but, nevertheless, shared a rather disputed 
simplification that perturbations in the flow thermophysical properties assume a single-harmonic 
type of circumferential variation. While such a simplification is, for all practiced purposes, applicable 
to the simple annular seal problem (Baskharone and Hensel, 1991b), it is unlikely that the same 
would hold true for secondary turbomachinery passages which, due to their complex geometry, give 
rise to flow separation and recirculation even under the unperturbed (centered-rotor) operation 
mode. 

The current study is an extension to a recent investigation (Baskharone and Hensel, 1991a), 
in which an unconventional perturbation approach to the general fluid-induced vibration problem 
was devised. The perturbation equations, under the new strategy, were deduced from the flow- 
governing equations in their discrete finite-element form, rather than the differential form as is 
traditionally the case. As a result, limitations on the circumferential perturbation pattern and/or 
the rotor-excitation degrees of freedom, were totally alleviated. The new perturbation model was 
utilized in computing the rotordynamic coefficients associated with cylindrical and conical rotor 


whirl [Baskharone and Hensel (1991b) & (1991c)], as well as a compound rotor excitation that is 
composed of the two whirling modes (Baskharone and Hensel, 1991d). 

Validation of the new perturbation model has so far been focused on verifying the results in 
a “macroscopic” sense, with the net integrated rotor forces, and the corresponding rotordynamic 
coefficients being the only variables under examination. In this paper, however, fine details of the 
perturbed rotor-to-housing flow field are rigorously verified. The intention here is to provide a 
better understanding of the source of rotordynamic forces which, in turn, should serve the design 
process itself. The current study was first and foremost motivated by the recent availability of 
detailed LDA flow measurements in a typical annular seal with a synchronously whirling rotor 
(Morrison et al., 1992). The outcome of this experimental study is utilized here for the purpose of 
comparison. 

The seal geometry and operating conditions in the current study are identical to those of 
Morrison et al. (1992). With a length of 37.3 mm, a rotor diameter of 164.1 mm, and a nominal 
clearance of 1.27 mm, this seal was tested at a rotor speed of 3600 rpm using water as the working 
medium. These variables gave rise to a Reynolds number (Re) of 24,000 and a Taylor number 
(Ta) of 6,600. At all times, the rotor whirl frequency was identical to the spinning speed, and the 
eccentricity ratio (^/c) was fixed at 50 % in the test rig. 

Computational Development. 

Centered-Rotor Flow Field. The perturbation procedure is initiated by computing 
this “zeroth-order” flow field, which is clearly axisymmetric. Details of the flow-governing equa- 
tions, turbulence closure, and finite-element formulation of this relatively simple problem were all 
discussed by Baskharone and Hensel (1991b), where two sample cases (including the seal under 
consideration here) were presented. Comparison of the through-flow and swirl velocity profiles in 


the rotor-to-housing passage with the experimental data previously reported by Morrison et al. 
(1991), was also presented, and the agreement between the two sets of data was encouraging. 

Perturbation Analysis. The centered-rotor flow solution, referenced above, was used 
in the current study as input to the perturbation analysis phase of the computational procedure. 
Devised by Baskharone and Hensel (1991a), the perturbation model here is based on what was 
generically termed the “virtually” deformable finite element concept, where the perturbed flow 
equations emerge from expansion of the finite-element equations in terms of the rotor eccentricity 
e. Since no circumferential pattern is pre-imposed on the perturbations of the flow properties, the 
problem at this point shifts to the fully three-dimensional type, as the rotor eccentricity destroys 
the flow axisymmetry (Fig. 2). Simultaneously, the flow problem is cast in a rotating-translating 
frame of reference, which is attached to the whirling rotor at all times (Fig. 2) in order to eliminate 
the time dependency of the flow field in this case. 

Grid Dependency. As would naturally be anticipated, the perturbed flow resolution 
is dependent on the grid resolution which, as the domain three-dimensionality prevails, translates 
into the number of computational stations in the circumferential direction (Fig. 2). Investigation of 
this dependency was carried out by varying this number and computing the fluid-exerted forces on 
the rotor each time. Results of this preliminary step are shown in Fig. 3, where the fluid force com- 
ponents F r ' and F e ’ are obtained by integrating the rotor-surface pressure forces, upon resolution 
in the radial and tangential directions, over all of the contributing finite element boundaries. 

The force trends in Fig. 3 suggest that the radial force is more sensitive to the tangential grid 
resolution by comparison. The Figure adso shows that changes in the force magnitudes becomes 
acceptably small for a tangential computational-station count that is in excess of 11. In the current 


study, this parameter was selected to be 13, as a compromise between a desirably high resolution of 
the flow properties in the circumferential direction, on one hand, and the CPU time consumption, 
on the other. 


Results and Discussion 

Examples of the “raw” results of the current study are shown in Fig. 4 for the seal middle 
cross section, as a representative axial location. These are contours of the non-dimensionalized 
perturbations in the velocity components; namely cV z ~ /V{ n , cV r /Vi„ and cV$ /Wj, where: 

Vi n is the seal-inlet average through-flow velocity 
u> is the rotor spinning speed 
Ti is the rotor radius 
c is the seal nominal clearance 

with the astrisk implying the operator In this, as well as all remaining Figures, the whirl fre- 
quency fl is identical to the rotor speed u>, giving rise to the case of a synchronous whirl investigated 
by Morrison et al. (1992). There is, however, no valid comparison between the results in Fig. 4 
and Morrison’s flow measurements, since the latter is a combination of the velocity perturbations 
and the unperturbed (centered-rotor) velocity field. 

Examination of Fig. 4 reveals that the maximum perturbation of the axial velocity component 
occurs in the vicinity of the rotor “pressure” side (where the clearance is smaller than the nominal 
value), and at a tangential location that lags the minimum-clearance position, by reference to the 
whirl direction. The radial velocity perturbations in this Figure are clearly smaller than those of the 
axial velocity, and attains its maximum value near the rotor pressure side as well. The difference, 
however, is that the the radial component peaks ahead of the minimum-clearance position in the 


whirl direction. As for the tangential velocity perturbation in Fig. 3, the peak value appears closer 
to the housing, and on the “suction” side relative to the rotor. 

Figure 5 shows contours of the non-dimensionalized pressure perturbations (p"c/ pV in 7 ) , as well 
as those of the combined (centered-rotor and eccentricity-related) pressure values in the annulus 
at the middle of the seal length. The combined pressure values in this Figure correspond to an 
eccentricity ratio (e/c) of 50% , and are non-dimensionalized using the seal-inlet static pressure 
(pin). As expected, the peak values, in each plot, occurs on the rotor pressure side, and not at the 
minimum clearance position as may be intuitively anticipated. Comparison of this observation, as 
well as those cited above in connection with the velocity perturbations, with the flow measurements 
by Morrison et al. (1992) was not possible, as the latter involved only the “net” velocity field in 
the whirling— rotor operation mode, with no pressure measurements whatsoever. 

A one— to— one comparison with Morrison’s experimental data are presented in Figs. 6 through 
10. In these Figures, the velocity component perturbations, corresponding to a rotor eccentricity 
ratio of 50%, are superimposed on the centered-rotor magnitudes. In assessing the numerical results 
in these Figures, it was both important as well as interesting to verify a fact that was solidly stated 
by Morrison et al., that the maximum through-flow velocity magnitude occurs on the pressure 
side of the rotor, over the early seed sections, but then moves towards the suction side as the flow 
progresses to the seal downstream sections. Comparison of the computed velocity contours in Figs. 
6 through 10 reveals that the numerical results are indeed in agreement with Morrison’s observation 

In interpreting the computed velocity components (e.g. Fig. 10), it may seem that the larger 
magnitude of the tangential velocity component “IV’ in the greater-clearance region is inconsistent 
with the principle of mass conservation. However, such observation would be valid should the 
problem at hand be of of such category as that where the axial mass flux is zero (e.g. the case 


where the flow is simply swirling beween two parallel endwalls, with no mass flux in the axial 
direction). The fact is that any valid control volume (in the current problem), to which the Vq- 
related mass flux is a contributor, has to have an axial extension between two cross-flow planes. 
In this case, the mass flux across these two planes has to also be taken into account in applying 
the continuity equation. With this in mind, and by reference to Figs. 8&10, let us (for instance) 
consider the left -hand-side half of the seal cross section, between the zj L — 0.77 and zjL = 0.99 
cross sections, as a control volume. As can be seen from the V z distribution over the left-hand-side 
halves of these two sections, neither the V z contours nor the average V z values are identical. The 
point can, therefore, be qualitatively made that the net mass flux due to what may appear as an 
invalid V$ distribution over the left -hand-side-halves of the annulus at the two cross sections, is 
balnced by that due to the The change in the V^-related flux between the two axial locations. 

It is important, in comparing the numerical and experimental sets of data, to point out that 
early sections of the seal (perhaps up to the mid-seal location) are those which no significant 
agreement would be anticipated. This is largely due to the flow admission losses which would 
prevail in the actual test rig, whereas the computational model treats this early seal segment 
as systematically, in the sense of turbulence closure and near-wall analysis, as does the entire 
seal. This is primarily why the velocity contours become more and more in agreement with the 
experimental data as the seal exit station is approached (e.g. Figs. 9 and 10). This, in particular, 
includes the tangential shift of the maximum axial-velocity position from a tangential location 
in the pressure (small clearance) region, gradually towards the suction side (Figs. 8 through 10) 
in a manner that is consistent with the flow measurements in these Figures. In fact, one would 
even go further, in reviewing the flow measurements in Figs. 6 through 10, indicating that the 
measurements themselves do not exhibit a characteristically similar axial-velocity pattern until the 
77% seal-length section is reached. 



An interesting phenomenon, which is confirmed both by the numerical results and the flow 
measurements, is that of flow separation and recirculation near the seal exit station. This, by 
reference to Figs. 9 and 10, takes place on the rotor suction side, and is more pronounced at the 
tangential location marked “s” in these two Figures. It is also in Figs. 9 and 10 that the numerical 
results yield shapes and magnitudes of the velocity contours which are very much consistent with 
the flow measurements. 

Another flow structure, which is equally interesting, involves the tangential velocity distribution 
and pertains, in particular, to the seal sections at 49% and 77% of the seal length (Figs. 7 and 
8). It is seen that there exists a region of reversed tangential velocity in each of these two sections 
near the housing at the minimum-clearance location. This implies local in-plane flow separation 
and recirculation in this region. Such flow behavior was apparently suspected (but not confirmed) 
by Morrison et al. (1992), who stated that: “if a tangential recirculation zone exists, it must be 
between the stator (referred to as the housing in this paper) and the first radial grid line.” 

Referring to the circumferential velocity contours in Figs. 8 through 10, in which the numerical 
and experimental results can be meaningfully compared, one can observe several characteristic 
differences between the two sets of data. Contrary to the numerical contours, the experimental V$ 
contours feature a local and abrupt fluctuation, midway between the rotor and the housing, at the 
greatest clearance position. As for the tangential location where the circumferential velocity peaks, 
one could see (by marching in the clockwise direction) that V$ peaks within the first quadrant in 
both sets of data. However, examination of the experimental contours in Figs. 8 through 10, does 
not reveal a fixed angular position where the peak value exists, although the angle (measured from 
the horizontal direction) of minimum V$ seems to fluctuate around an average value of 32°. This 
angle, by comparison, is nearly constant in the numerical set of data, and is approximately equal 


Referring to the numerical and experimental results in Figs. 6 through 10, it is important 
to point out that the uncertainties in the two sets of data are +6% and _+4.7%, respectively. 

Of these, the The experimental uncertainty was cited by Morrison et al. (1992). As for the 
numerical results, the uncertainty was deduced from the comparison by Baskharone and Hamed 
(1981) between the numerical results of a confined-flow problem and the closed-form exact solution 
of the same problem. The rationale here is that the finite-element formulation, the discretization 
unit, the matrix-inversion code and the numerical precision (in both computational models) are 
all identical to one another. 

Figure 11 shows a contour plot of the non— dimensionalized pressure perturbation; namely 
(p*c/ pVin) , and the resultant (centered-rotor and whirl-related) magnitudes of the non-dimension alized 
pressure ( pi Pin ) over the rotor surface. The pressure contours in this Figure are shown on the “un- 
wrapped” rotor surface which is obtained by “splitting” the surface along the lines labeled ti] 
and [a 2 -&2] m Fig- 1* Due to the lack of pressure measurements in Morrison’s study, no experi- 
mental verification of the results in Fig. 9 was possible. Examination of the pressure contours in 
this Figure reveals that the peak value occurs at a tangential location which lags the minimum- 
clearance position, by reference to the whirl direction. It is conceivable, and indeed likely, that the 
maximum pressure location is a function of the whirl direction (forward or backward) and/ or the 
the whirl frequency ratio (fi/w). For the synchronous whirl case under consideration, Fig. 9 shows 
that the rotor surface pressure peaks to a value of approximately 1.12 of the inlet pressure near the 
inlet station, and that with the streamwise loss in total pressure (due to friction) this magnitude 
gradually declines to approximately 0.93 of the inlet pressure at the exit station. Fig. 9 also shows 
that the maximum and minimum rotor-surface pressure locations hardly change their tangential 
location over the overwhelming majority of the seal length. 


Concluding Remarks 


Since the development of the so-called “virtually” deformable finite-element concept (Baskharone 
and Hensel, 1991a), the validation of this new approach never went beyond the traditional “macro- 
scopic” means of assessing the ultimate objective of the entire computational procedure; namely 
the final set of rotordynamic coefficients (Baskharone and Hensel, 1991b, c&d). Such a verification 
method is perhaps the least accurate, for it was always the global value of integrated rotor-surface 
pressure that was examined in reality. Recalling that the seal-clearance pressure field itself is 
among the least sensitive thermophysical properties, it is clear that the current comparative study 
is the most significant in appraising the mathematically rigorous computational model which centers 
around this categorically new finite-element-based concept. 

The numerical results corresponding to a synchronously whirling rotor of an annular seal are 
consistent with the flow measurements of Morrison et al. (1992). Agreement of the two sets of 
results is particularly favorable over the downstream half of the seal, which is where the real-rig 
flow admission losses have dissipated, and characteristically similar trends of the measurements 
themselves take place. Contrary to the common perception that perturbation analyses are only 
applicable to infinitesimally small rotor eccentricities, the favorable results of this study corresponds 
to an eccentricity that is as high as 50% of the seal nominal clearance. Successful validation of the 
computational model, under such conditions, underscores the model accuracy and establishes it as 
perhaps the most versatile perturbation approach, of its nature, in the general area of rotordynamics 

today. 
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Fig. 2 Distortion of the rotor-to- housing finite element assembly 
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Fig. 3 Relationship between the fluid-exerted forces and the grid resolution 
in the tangential direction 
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Fig. 8 Comparison of the velocity distribution at 77% of the seal length 
(eccentricity ratio = 50%) 







NUMERICAL RESULTS 


FLOW MEASUREMENTS 



Fig. 10 Comparison of the velocity distribution at 99% of the seal length 
(eccentricity ratio = 50%) 






